% $Id: all_au.tex 203 2009-04-20 01:20:13Z jsibert $
%
% Author: David Fournier
% Copyright (c) 2008 Regents of the University of California
%

%\tracingmacros=1
\documentclass[12pt]{book}
\usepackage{pictex}
\usepackage{makeidx}
%\usepackage{tocbibind}
%\usepackage{fancyhdr}
\input docmacpdf.tex
\input pdfcolor
%\usepackage{chappg}
%\usepackage{graphics}
\usepackage{hyperref}
\makeindex
\newcount\tm
\newcount\bm
\newcount\tom
\newcount\bom
\newtoks\myname
\newcounter{chapterpg}[chapter]
\pagestyle{empty}
\pagestyle{myheadings}
\makeatletter
\def\ps@yyy{
\def\@oddhead{}
\def\@evenhead{}
\def\@evenfoot{}
\def\@oddfoot{}
}

\def\ps@xxx{
\def\@oddhead{\message{Here DE}}
\def\@evenhead{\message{Here FE}}
\def\@oddfoot{
\plainline{ 
  \addtocounter{chapterpg}{1}
  \ifnum \number\c@chapter > 0
    \message{c@chapter = \number\c@chapter}
     \makebox[\textwidth]{
       {\ninerm admb-project.org}
      \hfill
     \arabic{chapter}-\arabic{chapterpg}\quad}
   \else
    \makebox[\textwidth]{\roman{chapterpg}}
   \fi
    \expandafter\write\newindex {%
      \thepage-\arabic{chapter}-\arabic{chapterpg}
     }
  %\makebox[\textwidth]{\arabic{page}}
 }} 
\def\@evenfoot{
\plainline{ 
  \addtocounter{chapterpg}{1}
  \ifnum \number\c@chapter > 0
    \message{c@chapter = \number\c@chapter}
       \makebox[\textwidth]{
       \arabic{chapter}-\arabic{chapterpg}\hfill
          {\ninerm {\slsmall Autodif} \ User's Manual}\quad}
   \else
    \makebox[\textwidth]{\roman{chapterpg}}
   \fi
    \expandafter\write\newindex {%
      \thepage-\arabic{chapter}-\arabic{chapterpg}
     }
  %\makebox[\textwidth]{\arabic{page}}
 }} 
}
\makeatother

\pagestyle{xxx}

\oddsidemargin=-.0truein
\evensidemargin=-.25truein
\textheight=8.5truein
\textwidth=6.5truein
\headheight=-0.0truein
\headsep=0pt
\topmargin=-.35truein
\topskip=-2.75truein
\begin{document}
     %\topmark=0
     %\botmark=0

\def\eqspace{\hskip 1000pt mi 990pt}

\def\endchapter{}
%\def\pageno{\thepage}
%\input ./docmacpdf.tex
\def\eqno{}
\def\thespace{ }
%\def\xfolio{\the\linknumber}
%\def\yfolio{\the\linknumber-\number\chapno-\number\pageno}
%\input pictex.tex
\advance\voffset by .5truein
\def\ADM{AD~Model~Builder}
%\def\chapno{1}
\font\rmmedium=cmr8 at 8truept
\font\rmfoot=cmr8 at 8truept 
\font\ttfoot=cmtt8 at 8truept 
%\font\ttsix=cmtt6 at 6truept 
%\font\stt=cmtt6 at 6truept 
\font\ttsix=cmtt8 at 6truept 
\font\stt=cmtt8 at 6truept 
\def\myeq#1{{\myeqwidth=\hsize \advance\myeqwidth by -.5in
   \hbox to \myeqwidth{\hfil$#1$\hfil}}}

\def\htmlnewfile{}
\def\htmlbeginignore{}
\def\htmlbegintex{}
\def\htmlendtex{}
\def\htmlendignore{}
%chapno=1
\thispagestyle{yyy}
\input adtitle.tex
\thispagestyle{xxx}
\tableofcontents
\newpage
\mychapter{Introduction}
\question{What is \AD}

\X{dvariable} \X{dvar\_vector}
\AD\ is an array language extension to \cplus\ which
allows the user to automatically calculate the partial derivatives
of a function of many independent variables by simply redeclaring the usual C floating point
types doubles and arrays of doubles to be the corresponding
\AD\ type, {\tt dvariable} and {\tt dvar\_vector}. 
The user can
write the code for the function using any valid \cplus\ syntax. 
\AD\ does not impose any restrictions on allowable code constructions. 
Using \AD\ the programmer can manipulate vectors and
matrices as single objects in complicated mathematical calculations,
at the same time obtaining derivatives of the resulting function with
respect to the independent variables.

When calculating derivatives there is no restriction on calling user supplied
functions which may themselves call other functions,
so long as the types have been redefined to be the corresponding \AD\ types.
Since \cplus\ is upwardly compatible with ANSI C,
existing numerical routines written in C can be easily 
modified to use \AD.

\question{Who should use \AD }
Anyone who would like to have the derivatives 
of a complicated function computed automatically
without any additional programming effort will find \AD\ useful.
 It is particularly useful for optimization problems involving
a differentiable function of many independent variables. Such
problems occur for example in nonlinear statistical modeling,
i.e. in nonlinear parameter estimation, in sensitivity analysis,
and in finding roots of systems of equations.

\question{Why is automatic differentiation better than 
estimating the derivatives by finite differences}
\X{automatic differentiation}

There are two reasons why the \AD\ method is better than finite differences
for calculating partial derivatives. 

\htmlbeginignore {\parindent=20truept\htmlendignore \XX{accuracy}{of derivative calculations}
\begin{enumerate}
\item With \AD\ the derivatives are calculated to the limit of machine
accuracy, that is to the same degree of accuracy as the function
itself. The derivatives are obtained as 
accurately as they would be if an analytical expression for the
derivatives were evaluated on the computer.
With finite differences the accuracy of the derivative
calculations seldom approaches that of the limit of the machine
accuracy for calculating the derivatives. 
Consequently, optimization schemes based on automatic
differentiation tend to be more stable and perform better than those
which use finite difference approximations.
 
\item To obtain estimates of the partial derivatives by finite difference
approximations the amount of calculation necessary
is proportional to the number of independent variables times the
amount of time needed to calculate the function itself. With \AD\ 
the amount of calculation is less than 5 times the amount of calculation 
needed to calculate the function itself (Griewank and Corliss 1991). 
In particular, the computing time does not depend on 
the number of independent variables.  This makes
\AD\ particularly attractive for calculating derivatives for
functions with many independent variables. 
\end{enumerate}
\par
\XX{speed}{of calculations}
\htmlbeginignore  }\htmlendignore 

\question{Is it necessary to learn \cplus before one can use \AD}

Any experienced C programmer can learn to use \AD\ in a few hours. While it
is necessary to have a \cplus\ compiler to use \AD,  only a small
number of extra concepts besides those used in the C language need
to be learned. This is due to the fact that \cplus\ is essentially
a superset of C, so that valid ANSI-compliant C code needs only a very small
amount of modification to be converted to valid \cplus\ code.

We feel that the best way to get a feel for how \AD\ is used is
to study a series of examples. These examples have been chosen to
illustrate the use of \AD\ on problems in different fields from
engineering to statistical analysis.

\question{How does automatic differentiation work}

Every computer program for calculating a function value, no
matter how complicated, can be broken down into a series of
binary and unary mathematical operations. Automatic differentiation
is simply the chain rule from elementary calculus applied to this 
series. The following simple example illustrates the chain rule.

%$$ \eqalign{ 
 \begin{eqnarray}
   u&=&x^2+y^2\cr
            v&=&x+3y\cr
            f&=&\sin(u)+\exp(v) 
 \end{eqnarray}
%}$$ 
\bestbreak
We want to calculate the derivatives of the function $f$ with
respect to the variables $x$ and $y$. By the chain rule
\def\pd#1#2{{\partial #1\over\partial #2}}

$$\pd fx=\pd fu\pd ux + \pd fv\pd vx\qquad\qquad
          \pd fy=\pd fu\pd uy + \pd fv\pd vy$$
Now
%$$\eqalign{
 \begin{eqnarray}
\pd fu=\cos(u)\qquad&\qquad\pd fv=&\exp(v)\cr
           \pd ux=2x\qquad\pd vx=1\qquad&\qquad&
           \pd uy=2y\qquad\pd vy=3
 \end{eqnarray}
%}$$
\noindent so that
$$\pd fx=\cos(u)*2x+\exp(v)*1\qquad\qquad\pd fy=\cos(u)*2y+\exp(v)*3$$
\noindent That is really all there is to it! Your program is
broken down into a series of simple operations like these
so that the derivatives can be calculated. The important
point is that no matter how complicated the
mathematical expressions are they can be thought of
as consisting of a long series of binary and unary operations.

\question{What is adjoint (precompiled) derivative code and 
  why is it so important}

While the reverse mode of automatic differentiation is
extremely efficient in terms of the number of calculations
required to compute derivatives, the method requires a
large amount of temporary storage. Each arithmetic operation 
used to calculate the original function generates 28--30 bytes
of temporary storage. Consider that to invert a 100 by 100
matrix requires $10^6$ operations. This would generate
28--30 megabytes of temporary storage requirements. Clearly
it would be impractical to carry out large calculations under 
these conditions. 
\X{precompiled derivative code}

The use of precompiled derivative code combines the best of two
worlds, so to speak. The small number of arithmetic operations
required by the reverse mode for calculating derivatives is
retained while at the same time  the amount of temporary storage
required is greatly reduced. The amount of storage required for inverting
a 100 by 100 matrix inverse is reduced to about 500K. Even better,
the amount of temporary storage requirement generated
by adding two matrices is about 32 bytes no matter how large
the matrices are.

Precompiled derivative code means that the actual derivative 
calculations are written in \cplus\ functions. Following 
object oriented design considerations the precompiled code is
completely ``encapsulated'' so that it is invisible to
the user. This means that it makes no difference to the user
whether a particular function has a precompiled derivative
code component written for it. At the user level everything looks the
same, except that performance is greatly enhanced.
The use of precompiled derivative code produces the ultimate in reuseable
code. For example, writing good code for the derivative of the 
the inverse of a matrix is not a trivial procedure, but once
it has been written any user of \AD\ can take advantage of it simply 
by taking the inverse of a matrix object with code like
\beginexample
M1=inv(M);  // M1 and M are dvar_matrix objects
\endexample

\AD\ combines array and matrix classes with the use of precompiled
derivative code for common array and matrix calculations and assignments
to produce a superior environment for efficient derivative calculations.

\section{A simple example of a derivative calculation}

While \AD\ has been supplied with all the computational power
described above, it has been designed so that at the user level
it appears like enhanced C code or perhaps as enhanced C code with
an adittional array language.  Before illustrating some of \AD's
array language facilities we will begin with a simple
C code example.

The following C code for the function {\tt fcomp}
calculates the sum of squares
$$\sum_{i=0}^{n-1} (x_i-1)^2 $$
This code has been written in C without using any \cplus\ or \AD\
extensions to show how C code can be modified to use \AD\ for calculating
derivatives.
\bigbreak
\beginexample
#include <math.h>
double fcomp(int n,double * x)
{
  double z;
  double tmp;
  int i;
  z=0;
  for (i=0;i<n;i++)
  {
    tmp=x[i]-1;
    z=z+tmp*tmp;
  }
  return(z);
}
\endexample
\bigbreak

\noindent{\tt fcomp} computes the sum of squares of n double precision floating
point numbers which are contained in an array of length {\tt n} pointed
to by {\tt x}. The expression {\tt x[i]} refers to the i'th
element of the array. (Note that the default construction for
C begins indexing an array with~0.) 

\section{Modifying fcomp to use \AD}
\X{dvariable}
The arithmetic 
 in the function {\tt fcomp} is carried out using
objects of type {\tt double}, that is, eight byte floating point
numbers. In order to calculate derivatives at the same time
as the function is being calculated, the class {\tt double}
is replaced with a new \AD\ class, the {\tt dvariable}. 
An object of type {\tt dvariable} ``appears'' to the user to
be an object of type {\tt double} in that it can be used in
arithmetic and mathematical expressions, and it has a numerical
value which can be displayed. The difference between an object of
type {\tt double} and an object of type {\tt dvariable} is that 
when arithmetic is done using a {\tt dvariable} object, information
necessary for the calculation of derivatives is automatically
generated as well. 

\X{dvar\_vector} 
In addition to objects of type {\tt dvariable}, we shall need
a class of ``vector'' object corresponding to {\tt x} which generates 
derivative information. In \AD\ this class is the {\tt dvar\_vector}
class. 


The following listing is the {\tt fcomp} function
modified to use automatic differentiation. 
\bigbreak
\listing{../../../examples/autodif/simple/fcomp_s.cpp}

\noindent There are five changes.

{\parindent=20truept
\begin{enumerate}
\item The line {\tt \#include <fvar.hpp>} has been included. The file
{\tt fvar.hpp}
contains information about \hbox{\AD}\ for the \cplus compiler.

\item The array of variables {\tt x} has been 
redefined to be of type {\tt dvar\_vector}. This \AD\ type is like
an array of doubles for which derivative information is collected. 
Besides providing for the calculation of derivatives, the use of the 
{\tt dvar\_vector} class for {\tt x} provides an additional improvement.
With the original C code there was no way of determining in the subroutine
{\tt fcomp} whether
enough memory had actually been allocated to the array~{\tt x}.
\AD\ provides optional index bounds checking for all container
classes such as {\tt dvar\_arrays}. This greatly speeds up 
program development.

\item The variables {\tt z} and {\tt tmp} have been declared to be
of type {\tt dvariable}.  The {\tt dvariable} type
is the fundamental type for automatic differentiation. It corresponds
to the C floating point type {\tt double}.

\item The function {\tt value} has been included
in the line {\tt return(value(z));}
The function  {\tt value()} returns the numerical part of an
{\tt dvariable}. This is necessary to convert the {\tt dvariable}
to a {\tt double} so that the value can be returned and used in
subsequent calculations (such as in a function minimizer) without
generating any more derivative calculations.

\item The range of the index of the {\tt for} loop has been changed
to begin at 1 and include n. While this change is not strictly necessary
fo this example,
\AD\ allows the user to declare arrays which have arbitrary valid
index ranges. This facility often contributes to the production 
of simnpler looking code.
\end{enumerate}

\par}

\section{Interfacing with the rest of the program}

The form of the interface of the function {\tt fcomp}
with the rest of the program will depend on the user's 
application. This section illustrates two applications.
The first example shows how to simply calculate the derivatives so that they can
be used for any other purpose whatsoever. 
The second example shows how to invoke a function minimizer, supplied with
\AD\ to
find the minimum of the function {\tt fcomp}. 
In any event it is necessary 
to identify the independent variables for \AD\ (these are the variables
for which the partial derivatives are calculated). The 
independent variables are identified by declaring them to be of type
{\tt independent\_variables}.\X{independent variables} 
The declaration {\tt independent\_variables x(1,nvar)} 
creates an array of 20 independent variables and identifies the 
independent variables for \AD. An object of type {\tt gradient\_structure}
must be declared and remain in scope until all the derivative calculations
have been finished.\XX{minimization}{of a function}\X{gradient\_structure} 
\beginexample
#include <math.h>
#include <fvar.hpp>
double fcomp(dvar_vector); // Function prototype declaration
void main()
{
  double f;
  int i;
  int nvar=20;
  independent_variables x(1,nvar);  // Identify the independent variables
  dvector g(1,nvar);  // Holds the vector of partial derivatives (the gradient) 
  gradient_structure gs;    // must declare this structure to manage derivative
                            // calculations
  f=fcomp(nvar,x);
  gradcalc(nvar,g);        // The derivatives are calculated
  cout <<" The gradient vector is\n"<<g<<"\n"; // Print out the derivatives
}                                              // on the screen
\endexample

\XX{constructor}{from {\tt independent\_variables} to {\tt dvar\_vector}}

Notice that although {\tt x} has been declared to be an object of
type\linebreak
{\tt independent\_variables} in the main routine, it is declared to
be an object of type {\tt dvar\_vector} in {\tt fcomp}. While this might
appear to be an error it is actually an important part of the
interface. It is necessary to identify both what the
independent variables are, and when one should begin to collect derivative
information for calculations involving them. When the vector
{\tt x} is passed to the function {\tt fcomp}, the compiler will
invoke a constructor to convert the type
{\tt independent\_variables} to the type {\tt dvar\_vector}.
At this time the steps necessary for
identifying the independent variables and initializing the
collection of information needed for calculating the derivatives 
 are carried out. With the Borland \cplus\ compiler
a warning message
\XX{Warning messages} {Temporary used for type xxx in function yyy}
\beginexample
Temporary used for type xxx in function yyy
\endexample
\noindent is printed
out notifying the user that a constructor has been invoked to perform
the necessary conversion. 
Such a mesage is usually benign and we tend to ignore them since we
have tried to design \AD\ so that unwanted type conversions do
not occur -- however it may indicate that some
unintended conversion has taken place. This is one of
the subtler pitfalls of \cplus\ programming. 
If you want to avoid such warning messages
it is possible to explicitly ``cast'' the object {\tt x} to the desired type
as
\XX{cast}{from {\tt independent\_variables} to {\tt dvar\_vector}}
\beginexample
  f=fcomp(nvar,dvar_vector(x));
\endexample
\X{BEGIN\_MINIMIZATION}\X{END\_MINIMIZATION}
\XX{macros}{BEGIN\_MINIMIZATION}\XX{macros}{END\_MINIMIZATION}
The next example shows how the above code is modified to 
invoke the function minimizer routine. The function to
be minimized and the function minimizer are embedded in a loop
so that it is not necessary for the user's function to call the
function minimizer or for the function minimizer to call the
user's function. This eliminates the need to communicate with
either routine by passing global variables.  Since the code which sets up 
the loop is independent of the user's function to be minimized we have
defined two macros, {\tt BEGIN\_MINIMIZATION} and {\tt END\_MINIMIZATION}
which will expand to produce the loop structure. The first macro, 
{\tt BEGIN\_MINIMIZATION} takes five arguments, the number of variables 
in the function to be minimized, the name of the vector of independent 
variables, the name of the dependent variable, the name of the 
vector of partial derivatives and the name of the {\tt fmm} control structure.
 It is the users
responsibility to declare these five objects.
The second macro, 
{\tt END\_MINIMIZATION} takes two arguments, the number of variables and
the name of the vector of partial derivatives. 
\bigbreak
\listing{../../../examples/autodif/simple/simple.cpp}

It is not necessary to use the macros. They have simply
been provided for convenience. If you prefer their expansion
can be written out as in the following example.
\beginexample
#include <math.h>
#include <fvar.hpp>
double fcomp(dvar_vector); // Function prototype declaration
void main()
{
  int nvar=20; // This is the number of independent variables
  independent_variables x(1,nvar); // these are the independent variables
  double f;  // This is the dependent variable
  dvector g(1,nvar);  // Holds the vector of partial derivatives (the gradient) 
  fmm fmc(nvar);      // Create structure to manage minimization
  // Expansion of BEGIN_MINIMIZATION(nvar,x,f,g,fmc) macro 
  gradient_structure gs;
  while (fmc.ireturn >= 0)   // Begin loop for minimization
  {
    fmc.fmin(f,x,g);         // Calls the minimization routine
    if (fmc.ireturn > 0)     // Loop for evaluating the function and
    {                        // derivatives
    // End of the Expansion of BEGIN_MINIMIZATION macro
      f=fcomp(nvar,x);
    // Expansion of END_MINIMIZATION(nvar,g) macro
      gradcalc(nvar,g); 
    }
  }
  // End of the Expansion of END_MINIMIZATION macro
  cout << " The minimizing values are\n" << x << "\n"; // Print out the answer
}
\endexample
The minimization routine is controlled by the object {\tt fmc} 
of type {\tt fmm}. Different aspects of the minimization, such as
the maximum number of function evaluations before stopping, convergence
criteria, and printing options can be adjusted by the user
through this object. See the chapter on the function minimizers for
more details. 
The flow through the loops is controlled by the {\tt integer fmc.ireturn}.

If you wish to use the conjugate gradient function minimizer instead 
of the quasi-Newton function minimizer you should declare the 
object {\tt fmc} to be of type {\tt fmmc} rather than type {\tt fmm}.
Nothing else needs to be changed.
\XX{function minimizer}{conjugate gradient}
\X{conjugate gradient function minimizer}

\section{Using \AD's built in vector and matrix calculations}

\XX{overloading}{of functions in C++} 
An important feature of \cplus\ is the ability of the language to
overload arithmetic operations, ie.  to extend their
definition to user defined types. For example the plus operation {\tt +}
can be extended to user defined string objects so that if
{\tt string1} and {\tt string2} are two string objects then
{\tt string1 + string2} produces a string object which contains
the concatenation of {\tt string1} and {\tt string2}. 
 \AD\ uses this feature to extend
the operations of addition, subtraction, multiplication, and division 
to various combinations of matrices, vectors, and numbers.
As a result the matrices and vectors can be employed as
objects in complex mathematical calculations while at the same time
derivatives are obtained. This greatly simplifies and speeds up
code development.
Addition and multiplication of vectors and matrices appear 
in many engineering and statistical applications.

\XX{bounds checking}{for container objects}  \X{safe arrays}
This example introduces two new types, the {\tt dvector} and the 
{\tt dmatrix}. These are the constant analogues of the {\tt dvar\_vector}
and the {\tt dvar\_matrix}. These are ``safe array''
objects which
contain their own size information and for which optional array 
bounds checking is provided.
The memory  for all these objects is
managed by \AD\ in a fashion which is transparent to the
user. \AD\ can combine these constant and variable objects automatically 
so that, for example, if {\tt A} is a {\tt dmatrix} and {\tt M}
 is a {\tt dvar\_matrix} the expression {\tt A*M} in the user code
will produce an {\tt dvar\_matrix} which is the product of {\tt A} and
{\tt M}. At the same time the array bounds on {\tt A} and {\tt M}
will be checked to ensure that the operation is legal for these two
objects. \X{dvector} \X{dvar\_vector} \X{dmatrix} \X{dvar\_matrix} 
 

To illustrate these vector and matrix operations we have taken the following 
example from Griewank (1988). 
The use of a cubic equation of state yields the 
Hemholtz energy
of a mixed fluid of unit volume at the absolute temperature $T$ as
$$ f(x)=RT\sum_{i=1}^n x_i \log\Big({x_i\over 1-b^Tx}\Big)
-{x^TAx\over \sqrt{8}b^Tx} 
   \log\bigg({1+\big(1+\sqrt{2}\big)b^Tx\over 
                                    {1+\big(1-\sqrt{2}\big)b^Tx}}\bigg)$$
where R is the universal gas constant and
$$0\le x,b \in {\hbox{\bf R}}^n\qquad,\qquad 
                      A=A^T\in {\hbox{\bf R}}^{n\times n}$$
During the simulation of an oil reservoir, this function and its
derivatives must be calculated at thousands of points in space and time.

In this example there are two new applications  of the overloaded multiplication operator ``{\tt *}''.
$b^Tx$ denotes the dot product of the  $n$ dimensional
vector~$b$ with the $n$~dimensional vector~$x$, 
$Ax$~denotes the  multiplication of the $n$~dimensional vector~$x$ by the
$n\times n$~matrix~$A$, and $x^T Ax$ is the result of the dot
product of the vector~$x$ with the vector~$Ax$. 
Here is the code for the function which calculates the 
Hemholtz energy using the overloaded operators.
\bigbreak
\listing{../../../examples/autodif/hemholtz/hem_func.cpp}

\X{overloaded operators}
All the vector and matrix multiplications are carried out using the
overloaded operators \hbox{\tt * + -} which enable the code to be developed 
extremely 
quickly.  Notice that the multiplication~{\tt *} can be carried out for
any combination of {\tt dvector}, {\tt dmatrix}, {\tt dvar\_vector}, and {\tt dvar\_matrix} so
long as the operation makes sense mathematically. \AD\ can combine the
constant types, {\tt double}, {\tt dmatrix}, and {\tt dvector}, with the variable types,
{\tt dvariable},
{\tt dvar\_vector}, and {\tt dvar\_matrix}, in a manner which is transparent to the
user.  

\goodbreak
Here is the the code which calls the function {\tt hemholtz\_energy}. 
\bigskip
\listing{../../../examples/autodif/hemholtz/hemholtz.cpp}

\section{The log-likelihood function for a multivariate normal distribution}
As a further example of to use \AD's vector and matrix
calculations consider the problem of estimating
the parameters for a multivariariate normal model
when the means and covariance matrix depend on a vector
of common parameters $\Theta$.
Let $Y$ be a vector of observations $(y_1, y_2,\ldots,y_n)$
which satifies the relationship
$$ Y = A(\Theta)X + \epsilon $$
where 
$X$ is a vector $(x_1, x_2,\ldots,x_m)$ of controls,
and $A(\Theta)$ is a $n\times m$ matrix depending on a 
parameter vector $\Theta=(\theta_1,\ldots,\theta_p)$,
and $\epsilon=(\epsilon_1, \epsilon_2,\ldots,\epsilon_n)$ 
is a random vector with a multivariate
normal distribution with mean vector $0$ and
covariance matrix $\Sigma(\Theta)$.

The log-likelihood function for the parameters $\Theta$
is given by
$$-0.5\ln\big(\hbox{\rm det}(\Sigma(\Theta))\big)-
   0.5(Y-A(\Theta)X)^T\Sigma(\Theta)^{-1}(Y-A(\Theta)X)$$
The maximum likelihood estimates for the parameters
are found by maximizing the
log-likelihood function with respect to $\Theta$. 
This is a nonlinear optimization problem. The best methods
for solving such a problem require the derivatives of
the log-likelihood function with respect to the parameter
vector $\Theta$. Using \AD's vector and
matrix classes, the code for calculating the log-likelihood
function and obtaining the derivatives can be written as follows. 
(Actually there are even more efficient methods of
coding these calculations with \AD.)
\XX{multivariate normal distribution}{log-likelihood function}
\bestbreak
{\openup\spread \openup\spread
\beginexample
dvariable log_likelihood(dvector x,dvar_vector Y,dvar_matrix Sigma,
  dvar_matrix A)
 // It is assumed that Sigma and A have already been calculated by
 // other routines in terms of the parameters Theta
{
  dvariable f;
  dvar_vector diff=Y-A*x;
  f=-0.5*log(det(Sigma))-.5*diff*inv(Sigma)*diff;
  return f;
}
\endexample}
When this code is executed, the information necessary for  calculating
the 
derivatives is automatically generated wiuthout further intervention
by the user. In addition optional array bounds checking is provided
to ensure that all vector and matrix objects have the right 
dimensions for the operations performed.

\question{What are ragged arrays}

Traditional programming languages such as FORTRAN which are often
used for numerically intensive applications only support
 ``flat'' data structures such as matrices. The natural
data structures for many applications simply do not come in this
simple form.  \AD\ has extended these simple forms to ragged matrices
and ragged three dimensional arrays.  A ragged matrix  is
a two dimensional array where the rows of the array can contain a 
different number of elements and have different valid index ranges.
A ragged three dimesnional array is made up of matrices of diferent 
sizes, and these matrices can themselves be ragged matrices.

Such ragged objects can be used to create many useful data structures.
For example the ``weights'' in a feed-forward neural network
can be best described by a ragged three dimensional array. This can be used to
provide an extremely compact description of the neural
network. See the chapter on neural networks for more details.

Since \cplus\ and \AD\ are extensible there
is in principle no limit to the complexity of the data structures
which can be developed. This means that you will never hit a ``dead end''
when developing an application.
\endchapter

\htmlnewfile 
%chapno=2
\mychapter{Getting Started}
\section{System Requirements -- PC implementations}
\X{System Requirements}
\X{hardware}
\AD\ is designed to operate on 80286, 80386 and 80486 based microcomputers
with a full 640 Kb main memory 
running under the MS--DOS opering system version 3.2 or higher.
A hard disk is
not required, however \AD\ may generate several Mb of gradient information
to be stored on an external device. A hard disk and a RAM disk are highly
recommended for greatest speed. The \AD\ libraries have been compiled with
either in line code for a numeric coprocessor or with code to emulate
a numeric coprocessor. A numeric coprocessor is also highly recommended
for greatest speed. 
The optimum configuration would include a numeric coprocessor plus a hard
disk and RAM disk with at least 2 Mb free memory for storing the temporary
files created by \AD.  (Large applications may require even more disk
memory.)
\AD\ appears to be compatible with most popular extended memory managers and
multitasking systems.

\X{libraries}
The \AD\ system consists of a header file {\tt fvar.hpp} and
two or more libraries. There is a safe library which provides bounds checking
and initialization of arrays. This library should be used for code
development. In addition there is an optimized libaray for faster execution
of production programs. Chekc the {\tt read.me} file on your \AD\
distribution disk for the names and number of libraries supplied with
your version of \AD.

\section{System Requirements -- other implementations}
\AD\ has been ported to a number of other platforms.
the necessary information for each specific platform should be
contained in a READ.ME file. 


\par{\parskip=0pt\parindent=6em
\par}

\section{Installation}

The simplest way to install the \AD\ system is to place the
header file in the directory where your compiler normally searches for
header files and the libraries in
the directory where your linker normally searches for
libraries.
If you choose this option, there is no need to specify additional
paths for the compiler and linker to search for header files and libraries.
Alternatively, you may wish to install the \AD\ system in its own
directory. In this case you will need to tell the compiler and linker
where to look for the header file and libraries.

An installation program is included on the distribution disk.
This program is invoked by typing {\tt a:install} at the DOS prompt
(assuming the distribution disk is in drive {\tt a:}).
You will then be asked for the drive of the source disk (in this case
{\tt a}) and the drive where you want to install \AD.
Next you will be asked for the directories in which to
install the \AD\ header file, libraries, and examples. 
If your \cplus\ compiler is located in {\tt c:\bs cc}, you might
want to specify {\tt c:} as the destination drive and
{\tt \bs cc\bs include} and {\tt \bs cc\bs lib} for the 
directories to 
install the include file and library files respectively.
If you press {\tt Enter} for
any of these directories the corresponding files will be installed in
a directory called {\tt \bs AUTODIF} on the destination drive.

{\bf Be sure to read the {\tt read.me} file for last minute changes and
compiler specific information.}

\section{Setting Up the Borland Compilers for \AD}
\X{Borland}
\XX{IDE}{Borland}
\X{memory model}
To compile \AD\ applications with the integrated development environment
(IDE), the following options must be selected within IDE:
\begin{description}
\item[---]Select {\bf O}ptions {\bf C}ompiler {\bf C}ode Generation.
Choose {\bf L}arge memory model. 
(Optional: Select Pre-c{\bf o}mpiled 
headers if you wish to speed up compilation.)
Select Mo{\bf r}e for Advanced Code generation. Select 8028{\bf 7} to
generate direct 80287 inline code or {\bf E}mulation to use numeric
coprocessor emulation. 
(Optional: Select {\bf D}ebug info in OBJs if you intend to use the
debugger.)
\item[---](Optional: If you intend to generate an overlay executable,
Select {\bf O}ptions {\bf C}ompiler {\bf E}ntry/Exit Code DOS {\bf o}verlay.)
\item[---]Select {\bf O}ptions {\bf D}irectories. If you have not stored
you \AD\ header file and library files in the default directories,
enter the path where you have stored the \AD\ header file ({\tt fvar.hpp})
under {\bf I}nclude Directories and the path where you have stored the
\AD\ libraries under {\bf L}ibrary Directories.
\item[---]Select {\bf P}roject {\bf O}pen to create a project file.
Include your source files and the \AD\ library that you wish to use.
\end{description}
\par

\X{make}
To compile \AD\ applications using a make file, the {\tt bcc} or {\tt bccx}
command-line compilers and {\tt tlink} or {\tt tlinkx} linkers must be used.

\XX{compiler options}{Borland}
\X{memory model}
You should invoke the following
compiler and linker options (see {\tt makefile} supplied with the examples):
\par{\parindent=4em\openup-\spread\openup -3pt
\begin{description}
\tracingmacros=1
\ttitem{-I}followed by the path of your \AD\ include file.
\tracingmacros=0
\ttitem{-ml}to specify the large memory model.
\ttitem{-L}followed by the path of your \AD\ libraries.
\ttitem{-c}to suppress linking after compilation.
\ttitem{-v}if you want to use the Borland debugger.
\ttitem{-f287} to use inline 80287 instructions.
\end{description}
\par}

\XX{linker options}{Borland}
You should invoke the {\tt /c} and {/v} linker options must be invoked
if your are using debugging.

\XX{linking}{Borland linkers}
The following files must be linked in the following order
(see {\tt makefile} supplied with the examples):
\par{\openup-\spread\openup -2pt
\begin{description}
\myitem{1.}{\tt c0l.obj} startup module for the large memory model.
\myitem{2.}Your source code {\tt .obj} files.
\myitem{3.}{\tt graphics.lib} if you have a graphic application.
\myitem{4.}{\tt ado7.lib} or {\tt ade.lib}
\myitem{5.}{\tt fp87.lib} for the inline floating-point library or 
{\tt emu.lib} for the emulation library.
\myitem{6.}{\tt mathl.lib} math library for large memory model. 
\myitem{7.}{\tt cl.lib} run-time library for large memory model.
\end{description}
\par}

\section{Setting Up the Zortech Compilers for \AD}
\X{Zortech}
\XX{IDE}{Zortech}
To compile \AD\ applications with the integrated development environment
(IDE), the following options must be selected within IDE:
\begin{description}
\myitem{---}Select {\bf C}ompile Com{\bf p}ile Options. 
Select {\bf M}emory Models L{\bf a}rge to select the large memory model.
Select Object {\bf C}ode Inline 808{\bf 7} to generate direct inline 8087 code.
\myitem{---}Append the directories where you have stored the \AD\ header
file and the \AD\ libraries to the environment strings {\tt INCLUDE} and
{\tt LIB} as discussed in the compiler documentation.
\end{description}

\XX{compiler options}{Zortech}
\X{memory model}
You should invoke the following
compiler and linker options (see {\tt makefile} supplied with the examples):
\par{\parindent=4em\openup-\spread\openup -2pt
\begin{description}
\ttitem{-I} followed by the path of your \AD\ include file.
\ttitem{-ml} to specify the large memory model.
\ttitem{-c} to suppress linking after compilation.
\ttitem{-g} if you want to use the Zortech debugger.
\ttitem{-f} if you want to generate in line 8087 code.
\end{description}
\par}

\section{Examples}
\X{examples}
Many of the examples in this manual are included on the distribution
disketted in a ``self-expanding'' archive program {\tt examples.exe}. When
this program is executed, it expands into  several files. 
The installation program will install the examples in the directory you
specify. A {\tt makefile} is included to simplify compilation of the
examples.
Read the file {\tt examples.doc} for a more complete description of the
examples.

\endchapter

\htmlnewfile 
%chapno=3
\mychapter{The \AD\ Classes}

\X{double}\X{float}\X{dvariable}\X{dvector}\X{dvar\_vector}
\X{constant objects} \X{container classes} \X{variable objects} 
\X{dmatrix}\X{dvar\_matrix}\X{d3array}\X{dvar3\_array}
\X{long int}\X{ivector}\X{lvector}\X{imatrix}\X{lmatrix}
\X{int}
The \AD\ floating point container classes have been designed to
facilitate the calculation of derivatives. As in differential calculus
they reflect
the main dichotomy between constant objects for which no derivative
calculations are required, and variable objects for which derivative
information is generated.  
Both constant and variable objects can be
further classified as single objects (a number), one dimensional arrays,
two dimensional arrays, and so on. The constant
and variable versions of each array type (container class) have
been designed to appear as similar as possible so that identical
code can be written for calculations involving either constant or
variable objects.

\vbox{\openup 3pt\tabskip=40pt plus40 pt minus 40pt
 \halign to \hsize{\hfil#&\hfil#&\hfil#\cr
 & Constant Type & Variable Type\cr
Single element&{\tt double}, {\tt float}&{\tt dvariable}\cr
1 dimensional array&{\tt dvector}&{\tt dvar\_vector}\cr
2 dimensional array&{\tt dmatrix}&{\tt dvar\_matrix}\cr 
3 dimensional array&{\tt d3array}&{\tt dvar3\_array}\cr} }%\vbox


\AD\ also has one and two dimensional container classes for integers.

\vbox{\openup 3pt\tabskip=40pt plus40 pt minus 40pt
 \halign to \hsize{\hfil#&\hfil#&\hfil#\cr
Single element&{\tt int}, {\tt long int}\cr
1 dimensional array&{\tt ivector}, {\tt lvector}\cr
2 dimensional array&{\tt imatrix}, {\tt lmatrix}\cr}}%\vbox 

The {\tt ivector} and {\tt lvector} 
classes are arrays of integers of type
{\tt int} and {\tt long int}.
The {\tt imatrix} and {\tt lmatrix} classes are two dimensional 
arrays of integers of type {\tt int} and {\tt long int}.

%\classdec{The {\tt dvariable} class}
\section{The dvariable class}

Objects of the {\tt dvariable} class are like numbers (variables).
When they are used in calculations, derivative information is generated. 
 The value of a {\tt dvariable} is
an object of type {\tt double}. All the standard C arithmetic and logical
functions have been defined for type {\tt dvariable}. 
A \hbox{dvariable {\tt z}} is created by the declaration
\beginexamplea
dvariable z;
\endexample

Since the {\tt dvariable} class is designed to appear to the user like
a {\tt double}  it is not necessary to discuss many of the operations which
can be performed. They appear to be simply the same as those performed
on an object of type {\tt double}.  There is however an important difference
between the behaviour of a {\tt dvariable} and that of a {\tt double}
under a copy operation. A copy occurs, for example, when an 
object is passed as an argument
to a function (so long as the function has not been declared to take
a reference). Suppose that the function {\tt func} has been declared
\X{copy operation} 
\beginexamplea
void func(dvariable)
\endexample
\noindent so that {\tt func} is a function which takes one argument of 
type {\tt dvariable} and a portion of the users code looks like
\beginexamplea
dvariable v;
// ...
func(v);
\endexample
\XX{constructor}{copy initializer}\X{dvariable(dvariable\&)}
When  the function {\tt func(v)} is called the copy initializer constructor
will be called first in order to construct a copy of {\tt v} and place it on 
the stack before passing control to
{\tt func}. For previous versions of \AD\ (before version 1.05) 
the copy initializer performed a 
shallow copy -- that is it
passed a pointer to the component of the dvariable containing the 
{\tt double}. As a result the {\tt dvariable} in the function {\tt func}
and the {\tt dvariable v} in the calling routine both point to the same
{\tt double} so that  changes to the value of the {\tt dvariable}
in {\tt func} will affect {\tt v} even though {\tt v} was passed by value.
To get better compatibility with the behaviour of C code (and the behaviour
of a {\tt double}) for version 1.05 or later the {\tt dvariable} copy initializer constructor has
been modified so that it produces a distinct copy. This could cause
a problem for old \AD\ code. If a {\tt dvariable} has been passed by
value to a function and its value modified in the body of the function
it used to be the case that the value of the {\tt dvariable} in the
calling function would be changed as well. This is no longer the case.
Care must be taken to modify old code where necessary. 
\XX{copy operation}{shallow}  

% 
% \XX{distinct copy}{using assignment operator = to create}  
% \XX{distinct copy}{of a dvariable}  
% This section applies only to versions of \AD\ with
% version number less than 1.05.
% For older versions of \AD\ (version number less than 1.05)
% when it is necessary to create a distinct copy of a {\tt dvariable} 
% the assignment operator~``{\tt =}'' should be used.
% To create a distinct copy of {\tt v} whose value is not affected by
% the function {\tt func} the preceding code could be modified to
% \beginexamplea
% dvariable v;
% dvariable vcopy; 
% vcopy=v; // make a distinct copy of v to pass to the function func
% // ...
% func(vcopy);
% \endexample
% \noindent
% {\bf It is important to realize that the following simpler code will not
% accomplish the objective of creating a distinct copy of {\tt v}.}
% \beginexamplea
% dvariable v;        // Not the right way to make a distinct copy
% dvariable vcopy=v;  // This will NOT make a distinct copy 
% // ...
% func(vcopy);
% \endexample
% \X{initialization}\XX{copy initializer constructor}{for a dvariable} 
% \noindent The problem is that while the code 
% {\tt dvariable vcopy=v;} appears to be an assignment, it is in fact 
% interpreted by the compiler as an initialization
% and will cause the copy initializer constructor to be invoked. 
% Thus a distinct copy of v is not created and the problem persists. 
% 

The question of distinct copies and the form of the copy initializer 
also occurs for container classes such as the
\AD\ vector and matrix classes. It is discussed below for the
{\tt dvar\_vector} class.

There is one other operation on {\tt dvariable}s which has no counterpart
for {\tt double}s, the \X{value function} value function. 
\beginexamplea
double value(dvariable&);
\endexample
\noindent This function returns a {\tt double} which is the value of 
the {\tt dvariable}. The main use of this function is to pass
the value of a some computation for which derivative information has been
calculated to some other part of the user's program so that the value
can be used without creating a derivative
calculation associated either with the passing of the value or its subsequent
use in other parts of the program.
It is not necessary to use {\tt value} for stream output 
of {\tt dvariable}s because the stream operators {\tt <<} and
{\tt >>} have been overloaded to accept these objects. 
The value function should be used sparingly. In particular it must not 
be used in the middle of calculations for which the derivatives are
desired. The following example will produce an error
\XX{stream operators}{<<\ \ >>} 
\XX{operators}{<<\ \ >>} 
\XX{stream input and output}{of dvariables}
\beginexamplea
dvariable x,y,z;
double a;
// ...             This code is incorrect
y=x*x;
a=value(y);
cout << " y = " << a << "\n";
z=sin(a); // This is where an error is introduced into the derivatives
cout << " z = " << z << "\n"; This will print out the value of z
\endexample
The following correct version of the code will produce exactly the same 
output and the derivative calculations will be correct 
\beginexamplea
dvariable x,y,z;
double a;
// ...             This code is correct
y=x*x;
cout << " y = " << y << "\n";
z=sin(y); // Use y instead of a at this point 
cout << " z = " << z << "\n"; This will print out the value of z
\endexample

It is possible to create an array {\tt u} of 10 {\tt dvariables} by 
using the declaration \XX{array}{of dvariables}
\beginexamplea
dvariable u[10]; // creates an array of 10 dvariables indexed from 0..9
\endexample
However, this is not the recommended construction. Instead the 
{\tt dvar\_vector} class should be used as in the declaration 
\beginexamplea
dvar\_vector u(0,9); // creates a dvar_vector of size 10 indexed from 0..9
\endexample
\noindent This construction generates more efficient code and produces
 a ``safe'' object with optional bounds
checking. The {\tt dvariable u[10]} construction may be useful for
converting pre-existing C code where the declaration {\tt double} or
{\tt float} is modified to {\tt dvariable}. 

\XX{float}{no corresponding \AD\ type}
There is at present no \AD\ type which corresponds to the C type {\tt float}. 
All \AD\ variable types contain a {\tt double}. 


\section{The dvector and dvar\_vector classes}

Since the {\tt dvector} and {\tt dvar\_vector} classes are intended to
appear identical to the user, with the exception that derivative
information is generated for the {\tt dvar\_vector} class, we shall
discuss the use of these classes simultaneously.

\X{safe arrays} \XX{maximum size}{of a dvar\_vector}
The {\tt dvar\_vector} class appears to the user as a vector (array)
of objects of type {\tt dvariable}. 
Bounds checking for array and matrix objects
can be implemented or not by linking the user's program with
the appropriate \AD\ library. The index for a {\tt dvector}
a {\tt dvar\_vector}
is of type {\tt int} (a signed integer) so that permissible array bounds
are -32,768 to 32,767 for 16 bit DOS versions where an object of type {\tt int} 
occupies two bytes.  The maximum number of elements which can be contained
in a {\tt dvar\_vector} is 8191 in 16 bit DOS versions of \AD.

\section{Creating dvectors and dvar\_vectors}
Several declarations can be used to create
a {\tt dvector} or {\tt dvar\_vector} object.
\beginexamplea
dvector v(1,n); // Creates a dvector with array bounds 1..n
dvector v=w;    // Creates a dvector and initializes it with w
dvector w("{1.2,4,-5.23e-1}");  // Creates a dvector with array bounds 1..3
                     // and initializes it with 1.2, 4, -5.23e-1.
dvector w("mystuff.dat"); //Creates a dvector whose contents are
                     // read in from the file mystuff.dat. The data in
                     // mystuff.dat must be invalid numeric format 
\endexample
The code {\tt dvector v=w;} creates a {\tt dvector v} 
by using  the ``copy initializer'' 
\penalty-10000 {\tt dvector(dvector\&)} 
constructor. This constructor performs a shallow copy so that
after this line of code is executed {\tt v} and {\tt w} will
both point to the same area of memory where the array is contained.
The same type of constructor is 
invoked when a {\tt dvector} or 
{\tt dvar\_vector} object is passed by value
to a function. For example suppose that a function {\tt users\_function}
has been defined where the function prototype is
\beginexamplea
void users_function(dvar\_vector v);
\endexample
\noindent and that the following statement appears somewhere in the user's code.
\beginexamplea
dvar\_vector w(1,n);
. . .
. . .
users_function(w);
. . .
\endexample
To pass the {\tt dvar\_vector w} to {\tt users\_function}, a copy of
{\tt w} is made and put on the stack. The copy initializer constructor
is called by the compiler to do this. 
The copy initializer uses a ``shallow'' copy to
create the copy of {\tt w} which is passed to the function.

A {\tt dvector} contains a pointer which points to the area of memory
where the elements of the {\tt dvector} are contained. With a shallow
copy the value of this pointer is simply passed to the new copy of
the {\tt dvector w}. As a result, both {\tt w} in the calling routine
and the copy of {\tt w} in the function {\tt users\_function} 
point to the same area of memory so that changes made to
this {\tt dvector} in {\tt users\_function} will change the values of
the corresponding entries of {\tt w} in the calling routine even
thought the {\tt dvector} was passed by value. This can create a
problem if the user wants to operate on the {\tt dvector} in 
{\tt users\_function} without changing the value of {\tt w} in the
calling routine. 

For ordinary C language objects,  
changes within a function to an argument passed by value to the
 function do not  affect the object passed.  This is not the case
for the \AD\ container classes {\tt dvariable, dvector, dvar\_vector, dmatrix,
or dvar\_matrix} or higher dimensional arrays. 
If you wish to pass a ``copy'' of an object to a function
so that the original object is not affected by changes to the
passed object you must explicitly make a distinct copy of the object
by using the assignment operator {\tt =} as in this example.
\XX{AVOID}{using a shallow copy when a distinct object is desired}.
\beginexamplea
dvar_vector v(1,n);  // want to pass a copy of v to the function func
dvar_vector w(1,n);  // This will be the distinct copy of v
w=v;                 // This creates the distinct copy
users_function(w);             // Can now call the function
\endexample
The user must avoid the simpler \cplus\ construction
of declaring an initialized object. Such a declaration
will cause the copy initializer to be invoked so that a shallow
copy will be performed.
\beginexamplea
dvar_vector v(1,n);  
dvar_vector w=v;     // this uses the copy initializer so that w is
                     // a shallow copy of v 
users_function(w);   // Changes to w in the function will cause changes 
                     // in v
\endexample

It is possible to create a dvector by explicitly naming the
data it will contain \XX{constructor}{dvector}
\beginexamplea
dvector w("{12,24,56 78 91.4 23.455e+2}"); 
\endexample
The data must be enclosed in double quotes and begin and end with
braces $\{$ $\}$. The numbers must be in valid numeric
format and be separated by blanks or commas. 
It is also possible to create a dvector using the contents of a file.
The declaration
\beginexamplea
dvector w("mystuff.dat");
\endexample
\noindent will create a {\tt dvector w} which contains the contents of
the file {\tt mystuff.dat}. It is not necessary to know how may numbers
are contained in the file to use this declaration. 
An error message will be generated if
the file does not exist, the file is empty, or there is non-numeric data in the
file. The data in the file must be separated by spaces. 
This constructor will create a {\tt dvector} whose minimum valid
index is 1.

It is also possible to create a {\tt dvector} or a {\tt dvar\_vector}
by reading in a column of data from a file. Suppose the file
{\tt ``work.prn''} contains the data
\medbreak
\vbox{
\tabskip=0pt plus 50pt
 \halign to \hsize{#\hfil&&\hfil #\cr 
Monday&  25.12& -41& 433.12&  10.11& 23\cr
Tuesday& 34.21&  356&  23&  23.1&  5\cr
Wednesday& 10&   3.13e-4& 43.11&  3.23&  4\cr 
Thursday& 12.1&   53.13& 453.1&  -5.13&  4.1\cr }}
\medbreak
The declaration
\beginexamplea
dvar_vector u("work.prn",3);
\endexample
\noindent will create a {\tt dvar\_vector} whose elements are the
third column of the file {\tt work.prn}. That is
\beginexamplea
u = (\ts-41,\ts356,\ts3.13e-4,\ts53.13\ts)
\endexample
\noindent The entries of the file
must be separated by blanks and all the entries in the column being
read must be in valid numeric format.
This constructor will create a {\tt dvector} whose minimum valid
index is 1.

% \section{Avoiding arithmetic statements in the initialization of
% {\tt dvariable} objects}
% {\bf This section only applies to version of \AD\ with
% version number less than 1.05} 
% 
% Including statements like
% \beginexample
% dvariable x=a*b;
% \endexample
% \noindent in your code can introduce errors. 
% This is due to the use of shallow copies in \AD\ combined with the
% the method used to allocate memory for objects the results containing
% %use of a circular queue of temporary objects  for containing
% intermediate calculations. If {\tt a} and {\tt b} are {\tt dvariables}
% then the operation {\tt a*b} is stored in a temporary variable.
% The shallow copy which is then used to initialize
% the variable {\tt x} will leave {\tt x} pointing to the area
% of memory allocated to the temporary variable. Later on this
% memory location may be overwritten leading to incorrect 
% results in the program. To avoid this type of error the proper
% construction is to split up this statement into two statements.
% \beginexample
% dvariable x;
% x=a*b;
% \endexample
% \XX{AVOID} {Initializing a {\tt dvariable} with an arithmetic expression}
% \noindent The statement {\tt x=a*b} performs a real copy so that
% {\tt x} does not share the same memory location as the temporary 
% variable.
% This restriction applies only to objects of type {\tt dvariable}.
% It is perfectly all right to use arithmetic statements in
% the initialization of variable array objects. For example
% \beginexample
% dvar_matrix M1=M*trans(M);
% \endexample
% \noindent where {\tt M} is a {\tt dvar\_matrix} or
% \beginexample
% dvar_vector v= 10*u;
% \endexample
% \noindent where {\tt u} is a {\tt dvar\_vector}.
% 

Since this version of \AD\ does not feature resizeable arrays,
the only way that you can place the producet of two 
{\tt dvar\_matrices}~{\tt M1} and~{\tt M2} in a {\tt dvar\_matrix}~{\tt M3}
 without knowing or 
calculating the size of the product is to use the construction
\beginexample
dvar_matrix M3=M1*M2;
\endexample
\noindent so that the copy initializer constructor will create 
an object with the correct size. Of course it is possible to calulate the 
size of the product and explicitly declare the {\tt dvar\_vector~M3}
by using the {\tt dvar\_matrix} access functions to obtain the
row and column dimensions of the matrices being multiplied  as in the code
\beginexample
dvar\_matrix M3(M1.rowmin(),M1.rowmax(),M2.colmin(),M2.colmax());
M3=M1*M2;
\endexample
\noindent but this solution lacks elegance.
\section{The dvector and dvar\_vector access functions}
Most of the components of the \AD\ classes are private so that 
they can not be accessed directly by the user. To allow the user 
to access the relevant elements of the classes,
public class member access functions are provided.
\beginexamplea
int indexmin();     //Returns the minimum allowable index  for an array
int indexmax();     //Returns the maximum allowable index  for an array
int size();         //Returns the size of an array
int shift(int min); //Changes the allowable indices to min and min+size-1
\endexample
These functions can be used to determine the size and valid index ranges of
a vector object which has been passed to a user's function. Thus
it is not necessary to explicitly pass such information when designing
a function which uses a vector object. The information is
already passed with the object itself.
If {\tt v} is a {\tt dvector}, the function {\tt v.indexmin()} will
return the minimum valid index for {\tt v}. The {\tt shift} function allows
the valid index range for a {\tt dvector} or {\tt dvar\_vector} to be
changed.
As an example of a function which
acts on a vector objects, here is a listing for a function which overloads
the multiplication operator {\tt *} so that {\tt v * w} will be 
the dot product (also sometimes called the scalar product)
of the {\tt dvectors v} and {\tt w}.
\beginexamplea
double operator * (dvector& t1, dvector& t2)
{
  // Check to see if vectors have the same valid index bounds
  if (t1.indexmin() != t2.indexmin() ||  t1.indexmax() != t2.indexmax())  
  {
    cerr << "Index bounds do not match in double operator * "
            "(dvector&,dvector&)\n";
    exit(1);  // Stop if vectors are not compatible
  }
  double tmp;
  tmp=0;
  for (int i=t1.indexmin(); i<=t1.indexmax(); i++)
  {
    tmp+=t1[i]*t2[i];
  }
  return(tmp);
}
\endexample
\XX{AVOID}{only creating a version of a function which takes a variable argument}

If you write your own functions it is a good idea to provide
both a constant and variable version at the time of writing.
This is especially true if you write a variable version of the
routine.
Suppose that the user has written a function 
\hbox{\tt void userfun(dvar\_matrix)}
which takes a \hbox{\tt dvar\_matrix} argument, but has neglected to
make a version \hbox{\tt void userfun(dmatrix)} which takes a 
{\tt dmatrix} argument. At some later time the user wants to
use the function {\tt userfun} with an argument of type
{\tt dmatrix} such as in the following code segment.
\beginexamplea
{
  dmatrix M(1,5,1,5);
  userfun(M); // M will be converted to a dvar_matrix
  // ...

}
\endexample
\noindent Since there is no version of {\tt userfun} which
accepts an argument of type {\tt dmatrix} the \cplus\ compiler
will convert {\tt M} to an object of type \hbox{\tt dvar\_matrix}
by invoking the \hbox{\tt dvar\_matrix(dmatrix\&)} constructor.
This is probably not what the user intends and will be disasterous
if an object of type {\tt gradient\_structure} is not in scope at the time. 

The corresponding variable version of the {\tt dot} product
for {\tt dvar\_vectors}
would look like \XX{dot product}{of vectors}
\beginexamplea
dvariable operator * (dvar_vector& t1, dvar_vector& t2)
{
  if (t1.indexmin() != t2.indexmin() ||  t1.indexmax() != t2.indexmax())  
  {
    cerr << "Index bounds do not match in dvariable operator * "
            "(dvar_vector&,dvar_vector&)\n";
    exit(1);
  }
  RETURN_ARRAYS_INCREMENT();
  dvariable tmp;
  tmp=0;
     
  for (int i=t1.indexmin(); i<=t1.indexmax(); i++)
  {
    tmp+=t1[i]*t2[i];
  }
  RETURN_ARRAYS_DECREMENT();
  return(tmp);
}
\endexample
\X{RETURN\_ARRAYS\_INCREMENT}\X{RETURN\_ARRAYS\_DECREMENT}
The appearance of the statements {\tt RETURN\_ARRAYS\_INCREMENT();}
and {\tt RETURN\_ARRAYS\_DECREMENT();} requires some explanation.
\AD\ employs a stack of queues to hold intermediate arithmetic
results for variable objects. To ensure that these queues are
not overflowed you should always increment the queue structure when
entering a function which returns a variable object and decrement
the queue structure when leaving the function.  If you fail to increment
the queue structure and a queue overflow occurs the program will
continue to operate, but the arithmetic will be corrupted.
If you fail to decrement the queue structure when leaving the function
you may eventually
increment past the end of the stack of queues  producing an error message
\beginexamplea
Overflow in RETURN_ARRAYS stack -- Increase NUM_RETURN_ARRAYS
There may be a RETURN_ARRAYS_INCREMENT()
which is not matched by a RETURN_ARRAYS_DECREMENT()
\endexample
Simply include the increment and decrement statements in every function
you create which returns a variable object and
everything will be fine. 
It is not necessary to include the increment and decrement statements
in a function which does not return a variable type because
such a function can not be used in an arithmetic statement.
{\bf You should never include the {\tt RETURN\_ARRAYS\_INCREMENT()} and
{\tt RETURN\_ARRAYS\_DECREMENT()} functions in a function which returns
a void or constant type such as {\tt double, dvector, or dmatrix}}.
The reason is that such a function might be called when there is
no object from the class {\tt gradient\_structure} in scope. 
Calling {\tt RETURN\_ARRAYS\_INCREMENT()} or 
{\tt RETURN\_ARRAYS\_DECREMENT()} when there is no object of
type gradient\_structure in scope could cause  
a serious program error.
\XX{AVOID} {Calling {\tt RETURN\_ARRAYS\_INCREMENT()}}  
\XX{AVOID} {Calling {\tt RETURN\_ARRAYS\_DECREMENT()}}  
\section{Accessing array elements}
\X{operators () []}
\XX{access functions}{array elements}
Both the {\tt ()} and the {\tt []} operators can be used to access elements 
of arrays. There is absolutely no difference between them. 
If {\tt v} is a {\tt dvector} the i'th
element of {\tt v} can be accessed by writing either {\tt v[i]}
or {\tt v(i)}.
Optional array bounds checking is provided for both operators.
Array bounds checking is implemented or turned off
by linking with the appropriate library. No changes in the user's code
are required. Linking with the safe libraries 
enables array 
bounds checking, while linking with the optimized libraries 
disables array bounds checking. 
\XX{libraries}{safe} 
\XX{libraries}{optimized} 
\section{Creating column vectors and row vectors}
It is important to recognize that there is a difference in \AD\ between
a vector object and a matrix object which has either 1 row or 1 column.
This means that a vector object is neither a row vector nor a column
vector, it is simply a vector object. We have adopted this
treatment of vector objects so that it is possible to
multiply vector objects such as {\tt u} and {\tt v}
simply by writing {\tt u*v}. If vector objects were the
same as matrices with 1 column, it would be necessary to
write {\tt trans(u)*v} to multiply two vectors. Here {\tt trans}
denotes the transpose of a {\tt matrix} object.

\X{column\_vector function}
\X{row\_vector function}\XX{constructor}{dmatrix}
Sometimes it is useful to be able to treat vector objects as
though they were matrices with one row or column. The functions
{\tt column\_vector} and {\tt row\_vector} enable this to
be done.  The function
\beginexample
dmatrix column_vector(dvector& v);
\endexample
\noindent will convert v into a matrix  with 1 column. The valid column
index is from 1 to 1. The valid row index bounds are the same 
as the valid index bounds for {\tt v}. The function
\beginexample
dmatrix row_vector(dvector& v);
\endexample
\noindent will convert v into a matrix  with 1 row. The valid row 
index is from 1 to 1. The valid column index bounds are the same 
as the valid index bounds for {\tt v}. 

As an example of how to use these function consider the problem
of computing the outer product of two vectors. The outer product
of two vectors $(u_i)$ and $(v_j)$ is the matrix $w_{ij}$
where $w_{ij}=u_iv_j$.  The statement
\beginexample
dmatrix w=column_matrix(v)*row_matrix(u);
\endexample
\noindent will compute the outer product of the {\tt dvectors}
{\tt u} and {\tt v}. Of course you can also use the
function {\tt outer\_prod} supplied with the AUTODIF libraries
to accomplish the same thing.
\beginexample
dmatrix w=outer_prod(u,v);
\endexample

\XX{operators}{<<} 
Another use for the {\tt column\_vector} function is to
print out a vector one element to a line. The overloaded
operator {\tt <<} will  write out a vector object all on
one line. This is useful if you want to put a vector object into a row
in a spreadsheet, for example, but it can be a problem for
some editors which can not read in such long lines in a file.
The code
\beginexample
cout << column_vector(v);
\endexample
\noindent will print out the vector v with one entry per line.
\section{The dmatrix and dvar\_matrix classes}
The {\tt dmatrix} and {\tt dvar\_matrix} classes
 are implemented
as arrays of {\tt dvector} and \hbox{\tt dvar\_vectors}. If {\tt M} is a
{\tt dmatrix} then {\tt M[j]} or {\tt M(j)} is a {\tt dvector}
whose elements are the {\tt j}th row of {\tt M}.
Since each row is a {\tt dvector} which contains its own size (shape)
information it is a simple matter to construct {\tt dmatrix} objects
which have a different numbers of rows in each column,
 so called ``ragged'' matrices. Of course these are not matrices in the
usual sense in which the word is used, 
but they can be useful container objects for
some applications so we have extended the  {\tt dmatrix} and
{\tt dvar\_matrix} classes
to include them.
As with the {\tt dvector and dvar\_vector} classes
{\tt dmatrix and dvar\_matrix} classes appear identical to
the user with the exception that derivative information
is collected for the {\tt dvar\_matrix class}.
An object of type {\tt dvar\_matrix} is created by using the declarations
\beginexamplea
dmatrix M(lbr,lur,lbc,luc); // Creates a (lur-lbr+1)$\times$(luc-lbc+1) matrix
dmatrix M1=M;               // M has already been defined. M1 is initialized 
                            // by M.
dmatrix M(lbr,lur,ivec_lb,ivec_ub); // Creates a ``ragged'' dmatrix with a 
                            // variable number of rows in each column
dmatrix M("filename");      // Creates a dmatrix whose elements are the 
                            // contents of the file ``filename''
dmatrix M( "{ 1.5,-2.0,3.0,3.5 }"
         "{ 2.0,14.0,5.1,2.2 }" );// creates and initializes a dmatrix
\endexample
The constructor {\tt dmatrix(lbr,lur,ivec\_lb,ivec\_ub)} 
creates a {\tt dmatrix} with {\tt lur-lbr+1} rows. The {\tt i}th row
is a {\tt dvector} of size {\tt ivec\_ub[i]-ivec\_lb[i]+1} and having
a range of valid indices running from {\tt ivec\_lb[i]} to
{\tt ivec\_ub[i]}. {\tt ivec\_lb} and {\tt ivec\_ub} are
instances of the class {\tt ivector}, that is vectors of
integers. They must have valid indices ranging from
{\tt lbr} to {\tt lur}. For example suppose we wish to make
a ragged {\tt dvar\_matrix} with 4 rows for which the valid row index bounds 
should be from 2 to 5. The declaration should be
\bestbreak
\beginexample 
dvar_matrix M(2,5,ivec_lb,ivec_ub);
\endexample
Suppose the desired valid row indices are 1 to 6 for row 1,
0 to 3 for row 2, 1 to 2 for row 4, and 3 to 9 for row 5.
The legal bounds on the {\tt ivector}s {\tt ivec\_lb} and {\tt ivec\_ub}
must be equal to the legal column bounds on {\tt M}. We can create such
{\tt ivector} objects with the declaration
\beginexample
ivector ivec_lb(2,5);
ivector ivec_ub(2,5);
\endexample
The correct values for {\tt ivec\_lb} and {\tt ivec\_ub} can
be inserted for example by using the ``fill'' function
\beginexample
ivec_lb.fill("{1,0,1,3}"); 
ivec_ub.fill("{6,3,2,9}"); 
\endexample
Alternatively the desired {\tt ivector} objects can be created with the
declarations
\beginexample
ivector ivec_lb("{1,0,1,3}"); 
ivector ivec_ub("{6,3,2,9}"); 
\endexample
These declarations will create {\tt ivector} objects whose legal bounds
go from 1 to 4. To change the valid index bounds to the desired values
from~2 to~5 the class member function {\tt shift()} can be used. 
\beginexample
ivec_lb.shift(2); // legal bounds will be from 2 to 5
ivec_ub.shift(2); // legal bounds will be from 2 to 5
\endexample

As is the case for vector objects, the copy initializer for matrix objects
performs a light copy. To create a distinct new matrix object you must
use the assignment operator {\tt =}.

The declaration
\beginexamplea
dmatrix M("filename");        
\endexample
\noindent will create a {\tt dmatrix M} whose elements are the contents 
of the file
{\tt filename}. 
All entries in the file
must be in valid numeric format. All entries must be separated by blanks
or commas.
Each line in the file becomes a line of the matrix. Ragged matrices are 
supported so that the files may have different number of entries in each 
line. For the DOS version of \AD\ no line may contain more than 6550
elements. 

To create a {\tt dmatrix} and initialize it with data a declaration of
the form
\beginexamplea
dmatrix W( "{ 1.5,-2.0,3.0,3.5 }"
           "{ 2.0,14.0,5.1,2.2 }" );
\endexample
\noindent can be used. This declaration will create a dmatrix with~2~rows
and~4~columns. The minimum valid index for the rows and columns is
set equal to~1.
\section{The dmatrix access functions}
\XX{dmatrix}{access functions}\XX{access functions}{for dmatrix and dvar\_matrix class}
The components of all the \AD\ classes are private which means that
they can not be accessed directly by the user. 
Public class member access functions are provided
to allow the user to access the relevant elements of the classes.
\beginexamplea
int rowmin();   //Returns the minimum allowable row index  for a matrix
int rowmax();   //Returns the maximum allowable row index  for a matrix
int colmin();   //Returns the minimum allowable column index  for a matrix
int colmax();   //Returns the maximum allowable column index  for a matrix
int rowsize();  //Returns the number of rows in an array
int colsize();  //Returns the number of columns in an array
void colshift(int&);  //Changes the range of valid indices for the columns
void rowshift(int&);  //Changes the range of valid indices for the rows 
\endexample
\X{rowmax()} \X{colmin()} \X{colmax()} \X{rowsize()}
\X{colsize()} \X{colshift(int&)} \X{rowshift(int&)}

\section{The three-dimensional arrays d3array and dvar3\_array}
\X{d3array and dvar3\_array classes}
Three dimensional arrays have been implemented in this version of
\AD\ mainly as containers for collections of two dimensional arrays. 
If {\tt u} is a {\tt d3array} then {\tt u[1]} or
{\tt u(1)} is a dmatrix (provided that 1 is a valid index for {\tt u}). 
We shall refer to the object
obtained by fixing the first index of a three dimensional
array as a ``slice''. An element of a three
dimensional array is determined by picking a slice, a row, and a column.
\XX{constructor}{d3array} \XX{slice}{of a d3array} 

An object of type {\tt dvar3\_array} is created by using the declarations
\beginexamplea
dvar_3darray M(ls,us,lr,ur,lc,uc); // Creates a (us-ls+1)*(ur-lr+1)*(uc-lc+1) 
                                   // dvar_3darray
dvar_3darray M1=M;                 // M has already been defined. M1 is 
                                   // initialized by M.
d3array(ls,us,ivec_lr,ivec_ur,ivec_lc,ivec_uc); // Creates a ragged d3array 
                            // with a variable size dmatrices in each slice.
\endexample
\section{The d3array access functions}\XX{access functions}{for a d3array}
\XX{d3array}{access functions}
\beginexamplea
int slicemin();   //Returns the minimum allowable slice index  for a d3array
int slicemax();   //Returns the maximum allowable slice index  for a d3array
int rowmin();   //Returns the minimum allowable row index  for a d3array
int rowmax();   //Returns the maximum allowable row index  for a d3array
int colmin();   //Returns the minimum allowable column index  for a d3array
int colmax();   //Returns the maximum allowable column index  for a d3array
int slicesize();  //Returns the number of slices in a d3array
int rowsize();  //Returns the number of rows in an d3array
int colsize();  //Returns the number of columns in an array
void sliceshift(int&);  //Changes the range of valid indices for the slices 
void colshift(int&);  //Changes the range of valid indices for the columns
void rowshift(int&);  //Changes the range of valid indices for the rows 
\endexample
\X{slicemin()} \X{slicemax()} \X{slicesize()}
\X{rowmin()} \X{rowmax()} \X{colmin()} \X{colmax()} \X{rowsize()}
\X{colsize()} \X{colshift(int&)} \X{rowshift(int&)} \X{slicehift(int&)}
\XX{conversion}{between between integer and floating point objects}
\section{Converting between the integer and floating point container objects}
It is often useful to be able to convert a floating point object into
an integer object or an integer object into a floating point object.
(See for instance the example on creating a bootstrap sample 
where a {\tt dvector} object
of random numbers is converted into an {\tt ivector} object of
integer indices).
It is possible to convert between vector objects as follows:
\beginexample
ivector <--> dvector
lvector <--> dvector
lvector <--> ivector
\endexample
These conversion may be made explicit as in the following code fragment
\beginexample
dvector u(1,10);
ivector iv(1,10);
// ...
iv = ivector(u); // iv will be filled with the integer parts of the components
                 // of u
\endexample
\noindent or the conversion may be ``implicit'', that is, supplied by the
compiler. \X{implicit conversion}\XX{conversion}{implicit}
\beginexample
dvector u(1,10);
// ...
ivector iv = u; // iv will be filled with the integer parts of the components
                // of u
\endexample
In the second example the {\tt ivector iv} 
will have minimum and maximum valid
indices determined by the minimum and maximum valid indices of the
{\tt dvector u}.

\endchapter


\htmlnewfile 
%chapno=4
\def\leftset#1{\quad #1 \hskip 1000pt minus 1fil}
\mychapter{Operations and Functions}
\section{The operators  + - * } \X{operators + - * }
\XX{vector}{component-wise sum of}\X{addition of vectors} 
\XX{optimizing performance}{using the best operators for a calculation}
Since \AD\ supplies a large number of operators and functions
there is often more than one way to carry out a given calculation.
While all equivalent methods will produce the same answer,
one method may produce more efficient code. For example,
the two lines of code
\beginexample
 A=A+B;
 //
 A+=B;
\endexample
\noindent will both add the object {\tt A}   to the object {\tt B}
and store the result in A. (This will work whether {\tt A} and {\tt B}
 are numbers,
vectors , or matrices.) However the code produced by {\tt A+=B}  
is more efficient in all cases. It is both faster and generates
less temporary data to be stored.

While the type of object produced by an arithmetic operation depends
on whether the arguments are constants ({\tt double,dvector or dmatrix})
or variables ({\tt dvariable, dvar\_vector, or dvar\_matrix})
the sort of operation defined by the operation does not.
The symbol {\tt +} between two vector objects will denote the component-wise
sum of these vectors and produce a vector object corresponding to their 
sum. This is true regardless of whether one or both objects are
constant or variable objects. We denote this by writing
\beginexample
vector_object = vector_object + vector_object // vector sum 
\endexample
We shall refer to the ``output'' object on the left side of the ``=''
sign by $z$ if it is a number, $z_i$ if it is a vector and $z_{ij}$
if it is a matrix in order to express the components of the
operation. We shall refer to the first input object (the one on the
right side of the ``='')
by $x$, $x_i$, and $x_{ij}$
and the second input object if there is one by
$y$, $y_i$, and $y_{ij}$
For any operation the property of being a variable is
dominant in that if either of the input objects
is a variable then the output object is a variable.
\beginexample
vector_object = vector_object + vector_object // vector sum 
\endexampledf
$$\leftset{z_i=x_i+y_i}$$
\beginexample
vector_object = number + vector_object // add number to vector 
\endexampledf
$$\leftset{z_i=x+y_i}$$
\beginexample
vector_object = vector_object + number // add number to vector 
\endexampledf
$$\leftset{z_i=x_i+y}$$
\beginexample
matrix_object = matrix_object + matrix_object // matrix sum 
\endexampledf
$$\leftset{z_{ij}=x_{ij}+y_{ij}}$$
\beginexample
matrix_object = matrix_object + number // matrix sum 
\endexampledf
$$\leftset{z_{ij}=x_{ij}+y}$$
\XX{vector}{component-wise difference of}\X{subtraction of vectors} 
\beginexample
vector_object = vector_object - vector_object // vector difference 
\endexampledf
$$\leftset{z_i=z_i-y_i}$$
\beginexample
vector_object = vector_object - number // subtract number from a vector 
\endexampledf
\XX{vector}{subtracting a number from a} 
$$\leftset{z_i=x_i-y}$$
\beginexample
matrix_object = matrix_object - number // matrix difference 
\endexampledf
$$\leftset{z_{ij}=x_{ij}-y}$$
\XX{matrix}{subtracting a number from a} 
\beginexample
matrix_object = matrix_object - matrix_object // matrix difference 
\endexampledf
$$\leftset{z_{ij}=x_{ij}-y_{ij}}$$
\XX{vector}{dot product of two}\X{dot product} 
\beginexample
number = vector_object * vector_object // vector dot product 
\endexampledf
$$\leftset{z=\sum_i x_iy_i}$$
\XX{vector}{multiplying a number by a}\X{scalar product} 
\beginexample
vector_object = number * vector_object // scalar product 
\endexampledf
$$\leftset{z_i=x*y_i}$$
\beginexample
vector_object = vector_object * number // scalar product 
\endexampledf
$$\leftset{z_i=x_i*y}$$
\XX{vector}{multiplying a matrix by a} 
\XX{matrix}{multiplying a vector times a} 
\beginexample
vector_object = vector_object * matrix_object // vector times matrix 
\endexampledf
$$\leftset{z_j=\sum_i x_i*y_{ij}}$$
\XX{vector}{multiplying a matrix times a} 
\XX{matrix}{multiplying a vector by a} 
\bestbreak
\beginexample
vector_object = matrix_object * vector_object // matrix times vector 
\endexampledf
$$\leftset{z_i=\sum_j x_{ij}*y_{j}}$$
\XX{matrix}{multiplication of two} 
\X{matrix multiplication} 
\beginexample
matrix_object = matrix_object * matrix_object // matrix times matrix 
\endexampledf
$$\leftset{z_{ij}=\sum_k x_{ik}*y_{kj}}$$
\XX{matrix}{multiplication by a number} 
\beginexample
matrix_object = number * matrix_object // scalar times matrix 
\endexampledf
$$\leftset{z_{ij}= x*y_{ij}}$$
An additional vector--matrix operation is the outer product of two vectors.
\XX{vector}{outer product of two} 
\XX{outer product}{of two vectors}
\beginexample
matrix_object = outer_prod(vector_object,vector_object) // outer product of two vectors
\endexampledf
$$\leftset{z_{ij}= x_i*y_j}$$
\section{The operators  += -= *= /=} \X{operators += -= *= /=}
These operators generate faster code than the corresponding operators defined
above since their use avoids the necessity for creating a 
temporary object to hold the result and the extra overhead associated with
assigning the values contained in this temporary object to the left hand side
of the expression. 
\beginexample
vector_object += vector_object // vector sum 
\endexampledf
$$\leftset{x_i+=y_i}$$
\beginexample
vector_object += number // add number to vector 
\endexampledf
$$\leftset{x_i+=d}$$
\beginexample
vector_object -= vector_object // vector difference 
\endexampledf
$$\leftset{x_i-=y_i}$$
\beginexample
vector_object -= number // subtract number from vector 
\endexampledf
$$\leftset{x_i-=d}$$
\beginexample
vector_object /= number // divide vector by number 
\endexampledf
$$\leftset{x_i/=d}$$
\beginexample
vector_object *= number // multiply vector by number  
\endexampledf
$$\leftset{x_i*=d}$$
\beginexample
matrix_object += matrix_object // matrix sum 
\endexampledf
$$\leftset{x_{ij}+=y_{ij}}$$
\beginexample
matrix_object += number // add number to matrix 
\endexampledf
$$\leftset{x_{ij}+=d}$$
\beginexample
matrix_object -= matrix_object // matrix subtraction 
\endexampledf
$$\leftset{x_{ij}-=y_{ij}}$$
\beginexample
matrix_object -= number // subtract number from matrix 
\endexampledf
$$\leftset{x_{ij}-=d}$$
\section{The operators { * /}} \X{* /}
\beginexample
vector_object = vector_object / number // divide a vector by a number
\endexampledf
\XX{vector}{dividing by a number} 
$$\leftset{z_i=x_i/y}$$
\beginexample
vector_object = number / vector_object // divide a vector by a number
\endexampledf
\XX{vector}{dividing  a number by} 
$$\leftset{z_i=x/y_i}$$
\beginexample
matrix_object = matrix_object / number // divide a matrix by a number
\endexampledf
\XX{matrix}{dividing by a number} 
$$\leftset{z_{ij}=x_{ij}/y}$$

\section{The concatenation of vector objects}
\XX{concatenation}{of vector objects} \XX{vector objects}{concatenation}
The {\tt \&} operator has been overloaded to produce the
concatenation of two vector objects.
\beginexample
vector_object=vector_object & vector_object
\endexample
\noindent If {\tt x}=$(x_r,\ldots,x_n)$
and {\tt y}=$(y_s,\ldots,y_m)$ then
{\tt x \& y}= $(x_r,\ldots,x_n,y_s,\ldots,y_m)$ and the minimum valid
index for {\tt x \& y} is equal to the minimum valid index for {\tt x}.

Since in the present version of \AD\ arrays are not resizeable
it is not possible to write something like
\X{resizeable arrays}
\beginexample
x = x & y;
\endexample
\noindent where {\tt x} and {\tt y} are {\tt dvectors}. Instead you
must do something like
\beginexample
dvector z = x & y;
\endexample
Future version of \AD\ will incorporate resizeable arrays.

\section{Element-wise operations}
There are several operations familiar to users of spreadsheets which do
not appear as often in classical mathematical 
calculations. For example spreadsheet users often
wish to multiply one column in a spreadsheet by the corresponding
elements of another column. Spread sheet users might find it much more
natural to define the product of matrices as an element-wise operation
such as
$$\leftset{z_{ij}= x_{ij}*y_{ij}}$$
The ``classical'' mathematical definition for the
matrix product has been assigned to the overloaded operator ``{\tt *}''
so that large mathematical formulas involving
vector and matrix operations can be written in a concise notation.
Typically, spreadsheet-type calculations are not so complicated
and do not suffer so much from being forced to adopt a
``function-style'' of notation.

Since addition and subtraction are already defined in an element-wise manner,
it is only necessary to define element-wise operations 
for multiplication and
division. We have name these functions {\tt elem\_prod} and {\tt elem\_div}.

\XX{vector}{element-wise product of}\XX{elem\_prod}{element-wise product}
\beginexample
vector_object = elem_prod(vector_object,vector_object) // element-wise multiply 
\endexampledf
$$\leftset{z_i= x_i*y_i}$$
\XX{vector}{element-wise division}\XX{elem\_div}{element-wise division}
\beginexample
vector_object = elem_div(vector_object,vector_object) // element-wise divide
\endexampledf
$$\leftset{z_i= x_i/y_i}$$
\XX{matrix}{element-wise product of}
\beginexample
matrix_object = elem_prod(matrix_object,matrix_object) // element-wise multiply 
\endexampledf
$$\leftset{z_{ij}= x_{ij}*y_{ij}}$$
\XX{matrix}{element-wise division}
\beginexample
matrix_object = elem_div(matrix_object,matrix_object) // element-wise divide
\endexampledf
$$\leftset{z_{ij}= x_{ij}/y_{ij}}$$
\section{The identity matrix function identity\_matrix}
\XX{matrix}{the identity matrix function}
\beginexample
matrix_object = identity_matrix(int min,int max)
\endexample
\noindent Creates a square identity matrix with minimum valid indices min and 
maximum valid index max.

\section{The operations  det inv norm norm2 min max sum}
\noindent The determinant of a matrix object (The matrix must be square, that 
is the number of row must equal the number of columns)
\XX{matrix}{the determinant of}\X{determinant}
\X{det}
\beginexample
matrix_object = det(matrix_object)
\endexample
\noindent The inverse of a matrix object (The matrix must be square, that 
is the number of row must equal the number of columns)
\XX{matrix}{the inverse of}\X{inverse} \X{\tt inv}
\beginexample
matrix_object = inv(matrix_object)
\endexample
\noindent The norm of a vector\_object
\XX{vector}{the norm of}
\XX{vector}{the norm squared of}\X{norm}
\XX{matrix}{the norm of}\X{norm}
\XX{matrix}{the norm squared of}\X{norm}
\beginexample
number = norm(vector_object)
\endexampledf
$$\leftset{z=\sqrt{\sum_{i} x_{i}^2} }$$
\noindent The norm squared of a vector\_object
\XX{matrix}{the norm of}\X{norm}
\beginexample
number = norm2(vector_object)
\endexampledf
$$\leftset{z=\sum_{i} x_{i}^2 }$$
The norm of a matrix\_object
\beginexample
number = norm(matrix_object)
\endexampledf
$$\leftset{z=\sqrt{\sum_{ij} x_{ij}^2} }$$
The norm squared of a matrix\_object
\beginexample
number = norm2(matrix_object)
\endexampledf
$$\leftset{z=\sum_{ij} x_{ij}^2 }$$
\goodbreak
\noindent The sum over the elements of a vector object
\XX{vector}{sum over the elements of}\XX{sum}{operation on a vector}
\beginexample
number = sum(vector_object)
\endexample
$$\leftset{z=\sum_i x_i }$$
\noindent The row sums of a matrix object
\XX{matrix}{rowsum of}\XX{rowsum}{operation on a matrix}
\beginexample
vector = rowsum(matrix_object)
\endexampledf
$$\leftset{z_i=\sum_j x_{ij} }$$
\noindent The column sums of a matrix object
\XX{matrix}{colsum of}\XX{colsum}{operation on a matrix}
\beginexample
vector = colsum(matrix_object)
\endexampledf
$$\leftset{z_j=\sum_i x_{ij} }$$
\noindent The minimum element of a vector object
\XX{vector}{minimum element of}\XX{min}{operation on a vector}
\beginexample
number = min(vector_object)
\endexample
\XX{vector}{maximum element of}\XX{max}{operation on a vector}
\noindent The maximum element of a vector object
\beginexample
number = max(vector_object)
\endexample
\section{Eigenvalues and eigenvectors of a symmetric matrix}
\XX{eigenvalues}{not differentiable}
\XX{eigenvectors}{not differentiable}
While we have included eigenvalue and eigenvector routines for
both constant and variable matrix objects you should be aware
that in general the eigenvectors and eigenvalues are not differentiable
functions of the variables determining the matrix.  

The eigenvalues of a symmetric matrix 
\XX{eigenvalues} {of a symmetric matrix}
\beginexample
vector_object = eigenvalues(matrix_object)
\endexample
\noindent are returned in a vector. It is the users responsibility to 
ensure that the matrix is actually symmetric. The routine symmetrizes 
the matrix so that the eigenvalues returned are actually those for
the symmetrized matrix.

The eigenvectors of a symmetric matrix 
\XX{eigenvectors} {of a symmetric matrix}
\beginexample
matrix_object = eigenvectors(matrix_object)
\endexample
\noindent are returned in a matrix. It is the users responsibility to 
ensure that the matrix is actually symmetric. The routine symmetrizes 
the matrix so that the eigenvectors returned are actually those for
the symmetrized matrix. The eigenvectors are located in the
columns of the matrix. The $i$'th eigenvalue returned by the
function {\tt eigenvalues} corresponds to the $i$'th eigenvector  
returned by the function {\tt eigenvector}.

\X{choleski decomposition of a symmetric matrix}
\XX{symmetric matrix}{choleski decomposition}
\section{The choleski decomposition of a positive definite symmetric matrix}
For a positive definite symmetric matrix {\tt S}, the
choleski decomposition of {\tt S} ia a lower triangular matrix {\tt T}
satifying the relationship {\tt S=T*trans(T)}. 
If {\tt S} is a (positive definite symmetric) matrix object and
{\tt T} is a matrix object, the line of code
\beginexample
T=choleski_decomp(S);
\endexample
\noindent will calculate the choleski decomposition of {\tt S} and put it
into {\tt T}.
\X{Solving a system of linear equations}
\X{solve function}
\XX{optimizing performance}{using the best operators for a calculation}
\section{Solving a system of linear equations}
If {\tt y} is a vector and {\tt M} is an invertible matrix then finding
a vector {\tt x} such that
\beginexample
x=inv(M)*y
\endexample
\noindent will be referred to as solving the system of linear equations
determined by {\tt y} and {\tt M}. Of course it is possible
to use the {\tt inv} function to accomplish this task but it is much
more efficient to use the {\tt solve} function.
\beginexample
vector x=solve(M,y); // x will satisfy x=inv(M)*y;
\endexample
It turns out that it is a simple matter to calculate the determinant
of the matrix {\tt M} at the same time as the system of linear
equations is solved, and since this is useful in multivariate
analysis we have also included a function which returns the 
determinant at the same time as the system of equations is solved.
To avoid floating point overflow or underflow when working with
large matrices the logarithm of the absolute value of the
determinant together with the sign of the determinant are returned 
The constant form of the solve function is
\beginexample
double ln_det;
double sign;
dvector x=solve(M,y,ln_det,sign);
\endexample
\noindent while the variable form is
\beginexample
dvariable ln_det;
dvariable sign;
dvar_vector x=solve(M,y,ln_det,sign);
\endexample
The solve function is useful for calculating the log-likelihood
function for a multivariate normal distribution.
Such a log-likelihood function involves a calculation similar to
\beginexample
  l = -.5*log(det(S)) -.5*y*inv(S)*y
\endexample
\noindent where Sis a matrix object and y is a vector object.
It is much more efficient to carry out this calculation using
the solve function. The following code illustrates the
calculations for variable objects. 
\XX{multivariate normal distribution}{calculation of 
  the log-likelihood function for}
\beginexample
   dvariable ln_det;
   dvariable sign;
   dvariable l;
   dvar_vector tmp=solve(M,y,ln_det,sign);
   l=-.5*ln_det-y*tmp;
\endexample

\section{Methods for filling arrays and matrices}
While it is always possible to fill vectors and matrices by
using loops and filling them element by element, this is tedious
and prone to error. To simplify this task a selection 
of methods for filling vectors and
matrices with random numbers or a specified sequence of numbers is available.
There are also methods for filling row and columns of matrices with
vectors. In this section the symbol {\tt vector} can refer to either
a {\tt dvector} or a {\tt dvar\_vector}.
while the symbol {\tt matrix} can refer to either
a {\tt dmatrix} or a {\tt dvar\_matrix}.

\X{filling arrays and matrices}
\beginexample
  void vector::fill("{m,n,...,}")
\endexample
\XX{vector}{fill}\XX{fill}{filling a vector} 
\noindent fills a vector with a sequence of the form
{\tt n, m, $\ldots$} The number of elements in the string must match the
size of the vector.

\XX{vector}{fill\_seqadd} \XX{fill\_seqadd}{filling a vector} 
\beginexample
  void vector::fill_seqadd(double& base, double& offset)
\endexample
\noindent fills a vector with a sequence of the form
{\tt base, base+offset, base+2*offset,$\ldots$}

For example if {\tt v} is a {\tt dvector} created by the statement
\beginexample
dvector v(0,4);
\endexample
\noindent then the statement
\beginexample
v.fill_seqadd(-1,.5);
\endexample
\noindent will fill {\tt v} with the numbers $(-1.0,-0.5,0.0,0.5,1.0)$.
\XX{matrix}{rowfill\_seqadd} \XX{rowfill\_seqadd}{filling a matrix} 
\beginexample
  void matrix::rowfill_seqadd(int& i,double& base, double& offset)
\endexample
\noindent fills row {\tt i} of a matrix with a sequence of the form
{\tt base, base+offset, base+2*offset,$\ldots$}

\XX{matrix}{colfill\_seqadd} \XX{colfill\_seqadd}{filling a matrix} 
\beginexample
  void matrix::colfill_seqadd(int& j,double& base, double& offset)
\endexample
\noindent fills column {\tt j} of a matrix with a sequence of the form
{\tt base, base+offset, base+2*offset,$\ldots$}

\XX{matrix}{colfill} \XX{colfill}{filling a matrix column with a vector} 
\XX{matrix}{rowfill}\XX{rowfill}{filling a matrix row with a vector} 
\beginexample
void matrix::colfill(int& j,vector&)
\endexample
\noindent fills the j'th column of a matrix with a vector 

\beginexample
void matrix::rowfill(int& i,vector&)
\endexample
\noindent fills the i'th row of a matrix with a vector 

In this section a uniformly distributed random number is assumed to have
a uniform distribution on $[0,1]$. A normally distributed random number
is assumed to have mean $0$ and variance $1$.
A binomially distributed random number is assumed to have a parameter $p$
where $1$ is returned with probability $p$ and $0$ is returned with
probability $1-p$.
A multinomially distributed random variable is assumed to have a vector
of parameters $P$ where $i$ is returned with probability $p_i$.
If the components of $P$ do not sum to $1$ the vector will be normalized
so that the components do sum to $1$.

\XX{vector}{fill\_randu} \XX{fill\_randu}{filling a vector with random numbers} 
\beginexample
void vector::fill_randu(long int& n)
\endexample
\noindent fills a vector with a sequence of uniformly distributed 
random numbers. The {\tt long int n} is a seed for the random number
generator. Changing {\tt n} will produce a different sequence of
random numbers.

\XX{matrix}{colfill\_randu} 
\XX{colfill\_randu}{filling a matrix with random numbers} 
\beginexample
void matrix::colfill_randu(int& j,long int& n)
\endexample
\noindent fills column {\tt j} of a matrix with a sequence of uniformly 
distributed random numbers
The {\tt long int n} is a seed for the random number
generator. Changing {\tt n} will produce a different sequence of
random numbers.

\XX{matrix}{rowfill\_randu} 
\XX{rowfill\_randu}{filling a matrix with random numbers} 
\beginexample
void matrix::rowfill_randu(int& i,long int& n)
\endexample
\noindent fills row {\tt i} of a matrix with a sequence of uniformly 
distributed random numbers

\XX{vector}{fill\_randbi} \XX{fill\_randbi}{filling a vector with random numbers} 
\beginexample
void vector::fill_randbi(long int& n, double& p)
\endexample
\noindent fills a vector with a sequence random numbers from
a binomial distribution.

\XX{vector}{fill\_randn} \XX{fill\_randn}{filling a vector with random numbers} 
\beginexample
void vector::fill_randn(long int& n)
\endexample
\noindent fills a vector with a sequence of normally distributed 
random numbers

\XX{matrix}{rowfill\_randn} 
\XX{rowfill\_randn}{filling a matrix with random numbers} 
\beginexample
void matrix::colfill_randn(int& j,long int& n)
\endexample
\noindent fills column {\tt j} of a matrix with a sequence of normally 
distributed random numbers

\beginexample
void matrix::rowfill_randn(int& i,long int& n)
\endexample
\noindent fills row {\tt i} of a matrix with a sequence of normally 
distributed random numbers

\XX{vector}{fill\_multinomial} \XX{fill\_multinomial}{filling a vector with random numbers} 
\beginexample
void vector::fill_multinomial(long int& n, dvector& p)
\endexample
\noindent fills a vector with a sequence random numbers from
a multinomial distribution. The parameter $p$ is a {\tt dvector} such that
{\tt p[i]} is the probability of returning $i$. The elements of {\tt p}
must sum to $1$.

\section{Methods for extracting from arrays and matrices}
\X{extracting data from arrays and matrices}
\XX{matrix} {extract\_column}\XX{extract\_column}{from a matrix}
\beginexample
vector extract_column(matrix& M,int& j)
\endexample
\noindent extracts a column from a matrix and puts it into a vector 

\beginexample
vector extract_row(matrix& M,int& i)
\endexample
\noindent extracts a row from a matrix and puts it into a vector.

\XX{extract\_diagonal}{from a matrix}
\beginexample
vector extract_diagonal(matrix& M)
\endexample
\noindent extracts the diagonal elements from a matrix and puts them into
 a vector.

\X{operator ()}
The function call operator {\tt (\ )} has been overloaded inv two ways to
provide for the extraction of a subvector. \XX{vector}{extracting a subvector}
\X{extracting a subvector}
\XX{matrix} {extract\_row}\XX{extract\_row}{from a matrix}
\XX{vector}{function call () to extract subvector} 
\beginexample
vector(ivector&)
\endexample
\noindent An {\tt ivector} object is
used to specify the elements of the vector to be chosen. If 
{\tt u} and {\tt v} are {\tt dvectors}  and {\tt i} is an {\tt ivector}
the construction
\beginexample
dvector u = v(i);
\endexample
\noindent will extract the members of v indexed by i and put them in the
{\tt dvector u}. The size of {\tt u}
is equal to the size of {\tt i}. The {\tt dvector u} will have
minimum valid index and maximum valid index equal to the minimum
valid index and maximum valid index of {\tt i}. 
The size of {\tt i} can be larger than the size of {\tt v} 
in which case some elements of {\tt v} must be repeated. The elements of
the {\tt ivector i} must lie in the valid index range for {\tt v}. 

If {\tt v} is a {\tt dvector} and {\tt i1} and {\tt i2} are two integers  
\beginexample
 u(i1,i2)
\endexample
\noindent is a {\tt dvector} which is a subvector of {\tt v} (provided of
course that {\tt i1} and {\tt i2} are valid indices for {\tt v}). Subvectors
can appear on both the left and right hand side of an assignment.
\beginexample
dvector  u(1,20);
dvector v(1,19);
v = 2.0;   // assigns the value 2 to all elements of v
u(1,19) = v; // assigns the value 2 to elements 1 through 19 of u
\endexample
\X{operator ++}
\X{operator --}
\XX{operator ++}{use with subvectors}
\XX{operator --}{use with subvectors}
\XX{operator ++}{for {\tt dvectors}}
\XX{operator --}{for {\tt dvectors}}
In the above example suppose that we wanted to assign the vector {\tt v}
to elements 2 through 20 of the vector {\tt u}. To do this we must first ensure that
they have the same valid index ranges. The operators {\tt ++} and {\tt --}
increment and decrement the index ranges by 1. The code fragment
\beginexample
dvector  u(1,20);
dvector v(1,19);
v = 2.0;   // assigns the value 2 to all elements of v
--u(2,20) = v; // assigns the value 2 to elements 2 through 20 of u
u(2,20) = ++v; // assigns the value 2 to elements 2 through 20 of u
\endexample
\noindent It is important to realize that from the point of view of the vector
{\tt v} both of the above assigments have the same effect. It will have
elements 2 through 20 set equal to 2. The difference is in the side
effects on the vector {\tt v}. The use of subvectors and increment and decrement
operations can be used to remove loops from the code. Note that
\XX{subvectors}{using to remove loops from code}
\beginexample
dvector x(1,n)
dvector y(1,n)
dvector z(1,n)
for (int i=2;i<=n;i++)
{
  x(i)=y(i-1)*z(i-1);
}
\endexample
\noindent can be written as
\beginexample
dvector x(1,n)
dvector y(1,n)
dvector z(1,n)
x(2,n)=++elem_prod(y(1,n-1),z(1,n-1));  // elem_prod is element-wise 
                                       // multiplication of vectors
\endexample
 
\section{Making a bootstrap sample}\XX{bootstrap sample}{code for}
As an example of how one can use some of these operations, consider the
problem of creating bootstrap samples from a set of data. 
We wish to create a random sample (with replacement) of the original data set.
Assume that the original data is contained in a {\tt dvector} object
named {\tt data}.
It is not necessary that the new samples be the same size as the original
sample so we shall allow the size of the bootstrap sample to be 
a variable. The code for the function {\tt bootstrap} follows.
\beginexample
dvector bootstrap(dvector& data, int& sample_size, long int seed)
{
  // data contains the data which is used to make the bootstrap sample
  // sample_size is the size of the desired bootstrap sample
  // seed is a seed for the random number generator 
  dvector tmp(1,sample_size);
  tmp.fill_randu(seed);  // 0<tmp(j)<1             for 1<j<sample_size
  tmp=tmp*data.size()+1; // 1<tmp(j)<data.size()+1 for 1<j<sample_size
  ivector iselect(tmp);  // Take the integer parts of tmp and put them
                         // into iselect
                         // 1<=iselect(j)<=data.size() for 1<=j<=sample_size
  tmp=data(iselect);     // Select the bootstrap sample and
                         // store it in tmp
  return(tmp);
} 
\endexample
The above  code could be employed in the user's code 
to make bootstrap samples by using something like the 
following code fragment 
\beginexample
// ...
dvector data(1,nobs)
// ... somehow the data is put into data
// ...
dvector boot_sample=bootstrap(data,200,1231); // Make a sample of size 200
                          // and create the dvector boot_sample to hold it
// ...
\endexample
\section{Sorting vectors and matrices}
\X{sorting}
\XX{dvector}{sorting a}
\XX{dmatrix}{sorting a}
\XX{sorting}{dvector}
\XX{sorting}{dmatrix}

While sorting is not strictly a part of methods for calculating the
derivatives of differentiable functions (it is a highly non-differentiable
operation) it is so useful for pre- and post-processing data that we
have included some functions for sorting {\tt dvector} and {\tt dmatrix}
objects. If {\tt v} is a {\tt dvector} the statement
\beginexample
dvector w=sort(v);
\endexample
\noindent will sort the elements of {\tt v} in ascending order and 
put them in the {\tt dvector} object {\tt w}. The minimum and
maximum valid indices of {\tt w} will be the same as those of {\tt v}.
If desired an index table for the sort can be constructed by passing
and {\tt ivector} along with the {\tt dvector}. This index tables can be
used to sort other vectors in the same order as the original vector
by using the {\tt ()} operator.
\X{operator ()}
\beginexample
dvector u={4,2,1};
dvector v={1,6,5}
ivector ind(1,3);
dvector w=sort(u,ind); // ind will contain an index table for the sort
// Now w=(1,2,4) and ind=(3,2,1)
dvector ww=v(ind);     // This is the use of the ( ) operator for subset
                       // selection. 
// Now ww=(5,6,1)
\endexample

The sort function for a {\tt dmatrix} object sort the columns of the
dmatrix into ascending order, using the column specified to do the
sorting. For example
\beginexample
dmatrix MM = sort(M,3);
\endexample
\noindent will put the sorted matrix into {\tt MM} and the third column of
{\tt MM} will be sorted in ascending order.

\section{Mathematical Functions}
 The following functions have been included in \AD\ by overloading the 
\cplus\ library functions
\beginexamplea
  sin cos tan asin atan acos sinh cosh tanh fabs (sfabs) exp log log10 sqrt pow 
\endexample
  \X{sin} \X{cos} \X{tan} \X{asin} \X{atan} \X{acos} \X{sinh} \X{cosh} \X{tanh} \X{fabs} \X{sfabs} \X{exp} \X{log} \X{log10} \X{sqrt} \X{pow} 
\noindent These functions can be used on numbers or vector\_objects
in the form 
\beginexample
number = function(number);
vector_object = function(vector_object);
\endexample
\noindent When operating on {\tt vector\_objects} the functions operate  
element by element, so that if {\tt y} is a {\tt dvector} whose
elements are $(y_1,\ldots,y_n)$ then {\tt exp(y)} is a 
{\tt dvector} whose elements are $(\exp(y_1),\ldots,\exp(y_n))$.

\XX{vector}{maximum element of}\XX{max}{operation on a vector}
\XX{vector}{minimum element of}\XX{min}{operation on a vector}
The functions {\tt min} and {\tt max} when applied to a {\tt vector\_object}
return a {\tt number} which is equal to the minimum or maximum element of the
{\tt vector\_object}

\section{Using the smoothed absolute value function sfabs}
The function {\tt sfabs} is a ``smoothed'' absolute value function
which agrees with the standard absolution function for all values
of its argument with absolute value greater than 0.001. For values
of its argument less than 0.001 in absolute value the function
uses cubic splines to produce a differentiable approximation to the
absolute value function. This function is intended to be used
when a differentiable function which approximates the absolute value
function is needed.

While the function {\tt sfabs} has been set to agree with the absolute
value function {\tt fabs} for all values greater than 0.001 
value has it
is simple matter to rescale the function to that it will
agree with the absolute value function over any desired range
by multiplying and dividing it by the desired scale factor.
\beginexample
double scale
double x
// ...
double v=scale*sfabs(x/scale); // This function will agree with fabs
                            // for all components of w with absolute value
                            // greater than scale*.001
\endexample
\bigbreak
$$\vbox{\input robgrph3.tex }$$
In figure \number\mychapno.1 the function {\tt sfabs} has been
``magnified'' by a scale factor of 1000 to show how it differs
from the ordinary absolute value function {\tt fabs}.
\bigbreak
\endchapter

\htmlnewfile 
%chapno=5
\mychapter{Advanced concepts}
This chapter is devoted to various matters which enable the
\AD\ user to create more advanced applications and to
obtain faster performance from the \AD\  code created.
Begining users  can skip or briefly skim the topics
contained here.
\XX{optimization considerations}{temporary files}
\section{Reducing the size of the temporary files}

A disadvantage of the reverse mode of automatic differentiation is
the large amount of temporary storage required for keeping
the intermediate data required for the derivative calculations.
 Every arithmetic operation
on a {\tt dvariable} object generates 28-30 bytes of temporary storage.
To see what this means consider that
to invert an $n\times n$ matrix requires about $n^3$ such operations,
so that the temporary storage required for inverting
a 100 by 100 matrix would be about 28-30 megabytes.

To reduce the requirements for temporary storage \AD\ employs
precompiled derivative code for common operation on variable
container classes. Depending on the operation this can greatly
reduce the temporary storage used. For the matrix inverse the
storage for a 100 by 100 matrix is about 500K, a reduction
of about 60 to 1. For a vector dot product the reduction in
temporary storage requirements is about 3.5 to 1. While
this is much less than the matrix inverse it can still be a considerable
saving in large problems.

To achieve the best performance you should
use \AD's vector and array operations on variable objects whenever
possible rather than writing your own code to manipulate elements of
these container classes. 
Consider the following code fragment which ``normalizes'' a
{\tt dvar\_vector} so that its components sum to {\tt 1}.
\beginexample
dvar_vector v(1,n);
// get some data into v somehow
// ...
// now normalize v so that its components sum to 1
dvariable sum=0.;
for (int i=1;i<=n;i++)
{
  sum+=v(i);
}
for (i=1;i<=n;i++)
{
  v(i)=v(i)/sum;
}
\endexample
\X{\tt sum}
\XX{optimizing performance}{using the best operators for a calculation}
The same normalization can be achieved by 
the following code which uses the \AD\ function {\tt sum} which
sums the components of a vector object as well as the divide operator  
for dividing a vector by a number. It will produce less temporary
derivative information.
\beginexample
dvar_vector v(1,n);
// get some data into v somehow
// ...
// now nomalize v so that its components sum to 1
v=v/sum(v);
\endexample
Even more efficient code is generated by employing the {\tt /=} 
operator which will divide a vector by a number.
\beginexample
v/=sum(v);
\endexample

\XX{optimization considerations}{Accessing container class elements}

\section{Creating efficient code when accessing arrays}
Consider the following code fragment for adding two matrix objects
\beginexample
dvar_matrix a(1,n,1,n);
dvar_matrix b(1,n,1,n);
dvar_matrix c(1,n,1,n);
// ...
for (int i=1;i<=n;i++)
{
  for (int j=1;j<=n;j++)
  {
    a(i,j)=b(i,j)+c(i,j);
  }
}
\endexample
A FORTRAN compiler which only has to deal with a small number of
array types can be expected to produce well optimized code
from the above source code. \cplus\ compilers are less likely
to be able to exploit that fact that the first index {\tt i} is constant 
in the inner loop. It is to be expected that in the future
the ability of \cplus\ compilers to optimize code will improve.
In the meantime code segments like the above can be
optimized by the user through the use of ``references''. 
If {\tt a} is a {\tt dvar\_matrix} then {\tt a(i)} 
(or {\tt a[i]}) is a {\tt dvar\_vector}. This fact can be used 
to reduce the redundant adressing operations on the first index
in the above example by defining reference objects for these
{\tt dvar\_vectors}. If you are not familiar with the
concept of references which are extremely useful in 
\cplus, consult your \cplus\ reference manual.
\beginexample
dvar_matrix a(1,n,1,n);
dvar_matrix b(1,n,1,n);
dvar_matrix c(1,n,1,n);
// ...
for (int i=1;i<=n;i++)
{
  dvar_vector& ai=a(i);
  dvar_vector& bi=b(i);
  dvar_vector& ci=c(i);
  for (int j=1;j<=n;j++)
  {
    ai(j)=bi(j)+ci(j);
  }
}
\endexample
\noindent Note however that the most efficient way to add two matrices
is to use the code
\beginexample
a=b+c;
\endexample
\bestbreak
Of course there is some overhead involved in establishing the references
so that this would not be worthwhile for very small values of {\tt n}.

\section{Ragged matrices and 3 dimensional arrays}
\XX{data structures}{ragged matrices}
\X{ragged matrices}
\XX{data structures}{ragged three dimensional arrays}
\X{ragged three dimensional arrays}
A big advantage of \cplus\ over traditional scientific 
programming languages such as FORTRAN is the increased
facility for defining and manipulating complex data structures.
Many real world situations can not be conveniently described in 
terms of the simple ``flat'' data structures represented by
vectors, matrices, and regular three dimensional arrays.

As an aid to creating more complex data structures \AD's
two and three dimensional
container classes have been extended to what we call
ragged arrays.  In short, a ragged matrix is a matrix whose
rows consist of vectors of different lengths and varying
valid index ranges, while a ragged three dimensional array
can be viewed as an array of matrices of different sizes
while the matrices which make up the array may or may not be
ragged themselves.

\section{Specifying ragged matrices}

To specify a ragged matrix requires two integers for the minimum
and maximum valid row indices, a vector of integers for the minimum
valid index of the vector forming each row of the matrix,
and a vector of integers for the maximum
valid index of the vector forming each row of the matrix. 
\beginexample
int min=0;
int max=4;
ivector minind(0,4);
ivector maxind(0,4);  
// get the desired values into the ivectors minind and maxind
// ...
dmatrix M(min,max,minind,maxind);
\endexample
Notice that the minimum and maximum valid indices for {\tt minind}
and {\tt maxind} must be the same as the minimum and maximum
valid row indices for the matrix {\tt M}.

Suppose that {\tt minind(3)=-1} and {\tt maxind(3)=5}. Then the 
fourth row of {\tt M} will be a vector with seven components
whose minimum valid index is {\tt -1} and whose maximum valid index
is {\tt 5}.

Often one desires to make a ragged matrix where the minimum valid index
for all the rows is the same. This can be accomplished by filling
the ivector {\tt minind} with the same number in each component.
For convenience a matrix constructor which takes an {\tt int}
for its third argument has been supplied. 
\beginexample
int min=0;
int max=4;
int minind=1;
ivector maxind(0,4);  
// get the desired values into the ivector maxind
dmatrix M(min,max,minind,maxind);
\endexample
\noindent This will make a ragged matrix {\tt M}. The minimum
valid indices  for vectors making up the rows of {\tt M}
will all be equal to {\tt 1}.

\section{Specifying ragged three dimensional arrays}
The range of constructors available for ragged three dimensional arrays
is considerably more complicated than for ragged matrices. 
The most general constructor for a ragged three dimensional array
has the following prototype.
\beginexample
d3_array(int sl,int su,ivector& rl,ivector& ru,imatrix& cl,imatrix& cu);
\endexample
The i'th slice of the {\tt d3\_array} is a ragged {\tt dmatrix} 
which is constructed using the data 
{\tt rl(i)\ ru(i)\ cl(i)\ cu(i)} as
described above. Notice that {\tt rl(i)\ ru(i)} are of type {\tt int} 
while {\tt cl(i)\ cu(i)} are objects of type {\tt ivector}
so that this all makes sense. This is the most ``ragged'' of
three dimensional objects. There is also a constructor with the following
prototype
\beginexample
d3_array(int sl,int su,ivector& rl,ivector& ru,ivector& cl,ivector& cu);
\endexample
\noindent The i'th slice of the {\tt d3\_array} will use the data 
{\tt rl(i)\ ru(i)\ cl(i)\ cu(i)} to make a {\tt dmatrix}.
Since these arguments are all of type {\tt int} each slice
of the ragged three dimensional array consists of a ``non-ragged''
or regular {\tt matrix} object.
An example of the use of a ragged three dimensional array for 
a data structure is given in the chapter on neural networks.
The complete list of constructors provided for specifying three dimensional
arrays is:
\beginexample
    d3_array( int& sl,  int& sh,  int& nrl,
       int& nrh, int& ncl,  int& nch);

    d3_array(int& sl,int& sh,ivector& nrl,ivector& nrh,
      imatrix& ncl,imatrix& nch);
 
    d3_array(int& sl,int& sh,ivector& nrl,ivector& nrh,
      int ncl,imatrix& nch);
 
    d3_array(int& sl,int& sh,ivector& nrl,ivector& nrh,
      ivector& ncl,ivector& nch);
 
    d3_array(int& sl,int& sh,int& nrl,ivector& nrh,
      int& ncl,ivector& nch);
 
    d3_array(int& sl,int& sh,int& nrl,ivector& nrh,
      int& ncl,int& nch);
 
    d3_array(int sl,int sh,int nrl,ivector& nrh,
      int ncl,imatrix& nch);
\endexample

\section{Complex data structures --- making order out of chaos}

As your model becomes more complex the number of \AD\
objects required will increase until the code again starts to become
unmanageable. One solution is to group these objects themselves
into structures or classes so that they can be manipulated as
single objects. 

Suppose that your model has 5 vector objects whose minimum valid
indices are~1 and whose maximum valid indices are all equal
to an integer {\tt n} whose value is known at run time. Denote
these vector objects by {\tt a1,a2,$\ldots$,a5}. Let the model
have 5 more vector objects whose minimum valid indices 
are~1 and whose maximum valid indices are {\tt 2n}.  
Denote these vector objects by {\tt b1,b2,$\ldots$,b5}.
Finally suppose the model has a (regular) matrix object whose valid
row indices go from {\tt 1} to {\tt n} and whose valid column
indices go from {\tt 0} to {\tt 2n}. Denote this matrix by {\tt M}.
We can encapsulate all these objects in a class as follows.
\beginexample
// class definition for constant_model_parameters
class constant_model_parameters
{
public:
  int size;
  dvector a1; 
  dvector a2; 
  dvector a3; 
  dvector a4; 
  dvector a5; 
  dvector b1; 
  dvector b2; 
  dvector b3; 
  dvector b4; 
  dvector b5; 
  dmatrix M;
  constant_model_parameters(int n); // prototype for constructor 
};
\endexample
At the point in your program where you want to ``create'' this object
you will need to invoke the constructor described below. The invocation
is simply a declaration of the object with its name such as:
\beginexample
int n=10;
// ...
constant_model_parameters cmp(n); //creates the object with name "cmp"
\endexample
The main part of the implementation of the class is in the definition
of the constructor.  
\beginexample
 constant_model_parameters::constant_model_parameters(int n) :
   a1(1,n),a2(1,n),a3(1,n),a4(1,n),a5(1,n), 
   b1(1,2*n),b2(1,2*n),b3(1,2*n),b4(1,2*n),b5(1,2*n), 
   M(1,n,0,2*n) {size=n;}
\endexample
\noindent If you are not familiar with the form of the constructor above
you should consult a \cplus\ reference book. The idea is that the
constructors for all the classes which are members of
the class {\tt constant\_model\_parameters}
are invoked after the {\tt :} and before the body of the constructor.

Note that while for the ``safe'' library all
container classes are initialized to {\tt 0}, for the optimized
library they are not. If you want to initialize some of these
objects, such as {\tt a1} and {\tt M} this should be done
in the body of the constructor.
\beginexample
 constant_model_parameters::constant_model_parameters(int n) :
   a1(1,n),a2(1,n),a3(1,n),a4(1,n),a5(1,n), 
   b1(1,2*n),b2(1,2*n),b3(1,2*n),b4(1,2*n),b5(1,2*n), 
   M(1,n,0,2*n) 
 {
   size=n; 
   a1.initialize();  // set components of a1 and M to 0
   M.initialize();
 }
\endexample
All the objects in the class {\tt constant\_model\_parameters}
are as the name implies, constants. This class could be useful for reading in 
the parameters from a file, for example, and otherwise manipulating 
them. At some point presumably the user wants the parameters
transformed into variable objects so that derivative information
can be generated. This can be done by defining a corresponding 
variable class.

\beginexample
// class definition for variable_model_parameters
class variable_model_parameters
{
public:
  dvar_vector a1; 
  dvar_vector a2; 
  dvar_vector a3; 
  dvar_vector a4; 
  dvar_vector a5; 
  dvar_vector b1; 
  dvar_vector b2; 
  dvar_vector b3; 
  dvar_vector b4; 
  dvar_vector b5; 
  dvar_matrix M;
  // This constructor transforms the constant parameters into 
  // variable parameters
  variable_model_parameters(constant_model_parameters&); // prototype 
                                                     // for constructor 
};
\endexample
The constructor  for the class {\tt variable\_model\_parameters}
takes the constant objects in the class {\tt constant\_model\_parameters} 
and converts them to variable objects.
\beginexample
 variable_model_parameters::variable_model_parameters
   (constant_model_parameters& cms) :
   a1(cms.a1),a2(cms.a2),a3(cms.a3),a4(cms.a4),a5(cms.a5), 
   b1(cms.b1),b2(cms.b2),b3(cms.b3),b4(cms.b4),b5(cms.b5), 
   M(cms.M) {size=cms.size;}
\endexample
After all this has been done the constant and variable classes could be 
used something like
\beginexample
void fcomp(variable_model_parameters&); // function prototype for fcomp
int n=10;
constant_model_parameters cms(n);
// get some data into cms
//
// call a function which takes an object of type variable_model_parameters
fcomp(cms); // cms will be converted into an object of type
            // variable_model_parameters
\endexample
\noindent Since in this example, the only constructor which has been
defined for an object of type 
{\tt variable\_model\_parameters} takes an argument of type
{\tt constant\_model\_parameters} the only way to create the variable object
will be to first create the constant object. This may or may not be
what is required. Of course a constructor  can
be defined which will create a variable object  from the {\tt int n}.

\section{Default constructors for complex data structures}
You should be aware that the \cplus\ compiler attempts to 
create default constructors for complex classes if they are not provided by the
user. To do this it employs the constructors provided for the members of 
the complex class. For examples if you employ the lines of code

\beginexample
constant_model_parameters cmp(10);
constant_model_parameters cmp1=cmp;
\endexample

\noindent the compiler will construct a default copy constructor
for the class 
\linebreak {\tt constant\_model\_parameters} (assuming that you,
the user, have not already defined one).
In order to create the default copy constructor the compiler will invoke the 
copy constructors provided by \AD\ for {\tt dvector} and
{\tt dmatrix} classes. Since these constrtuctors perform
shallow copies, the {\tt dvector} and {\tt dmatrix} 
objects in  {\tt cmp1} will not be distinct from the corresponding
objects {\tt cmp}. If you want to define a copy constructor for
the class {\tt constant\_model\_parameters} which will create
a distinct object it could be done as follows:
\beginexample
 constant_model_parameters::
   constant_model_parameters(constant_model_parameters cms) :
   a1(1,cms.size),a2(1,cms.size),a3(1,cms.size),a4(1,cms.size),
   a5(1,cms.size),b1(1,2*cms.n),b2(1,2*cms.n),b3(1,2*cms.n),
   b4(1,2*cms.n),b5(1,2*cms.n), 
   M(1,cms.size,0,2*cms.size) 
 {
   size=cms.size;
   a1=cms.a1; 
   a2=cms.a2; 
   a3=cms.a3; 
   a4=cms.a4; 
   a5=cms.a5; 
   b1=cms.b1; 
   b2=cms.b2; 
   b3=cms.b3; 
   b4=cms.b4; 
   b5=cms.b5; 
   M=cms.M;
 }
\endexample



\endchapter
\htmlnewfile 
%chapno=6
\mychapter{The \AD\ libraries}

\XX{AUTDIF}{libraries}\XX{libraries}{safe}\XX{libraries}{optimized}
\AD\ is supplied with safe and optimized versions
 of its object module libraries.
The safe version is intended to be used for code development.
It provides for bounds checking of arrays and initialization
of all declared objects. This initialization sets the
value of all array members to 0 and initializes the derivative
structure associated with each variable type object.
The initialization of the derivative structures can
involve a fair amount of overhead so that after a project
has been debugged the user may wish to use the optimized version
of the library. The optimized version does not provide array
bounds checking or initialization of objects. 
It is recommended that the user initializes all objects which require
initialization during
the development stage of the program even though this is not
necessary when using the safe library so that the project code can
be easily linked with the optimized library later if desired.
If an application runs successfully using the safe library, but does
not run properly using the optimized library, there is probably
an unitialized object being used somewhere in the code. 

The names and number of libaries supplied with \AD\ vary from
compiler to compiler, so check the {\tt read.me} file on
the distribution disk for the names of the libraries supplied with
your version of \AD.

\X{initializing AUTODIF classes}
The class member functions {\tt initialize} are provided for initializing all
appropriate \AD\ classes. For example to declare and 
initialize a dvar\_matrix M you could use the following code fragment.
\beginexample
dvar_matrix M(1,20,1,30) // Creates a 20 by 30 matrix
M.initialize();          // Sets all the elements of M to 0 and initializes
                         // the derivative structures for M
\endexample
An example of code where an object requires initialization is 
\beginexample
{
  // . . .
  dvariable u; //u is declared but not initialized
  for (int i=1;i<=n;i++)
  {
    u+=x[i];  // x is a vector object which was previously declared
  }            // This is an error for the optimized library
 }
\endexample
\bestbreak
In this example the {\tt dvariable u} is declared, but not initialized.
The use of the unitialized {\tt u} in the expression {\tt u+=x[i]}
is an error if the code is linked with the optimized library.
It will perform fine if the code is linked with the safe library.
To use the optimized library this code could be rewritten as
\beginexample
{
  // . . .
  dvariable u=0; //u is declared and initialized
  for (int i=1;i<=n;i++)
  {
    u+=x[i];  // x is a vector object which was previously declared
  }
}
\endexample
\noindent The code could also be written as
\beginexample
{
  // . . .
  dvariable u; //u is declared but not initialized
  u=0;         //u is initialized
  for (int i=1;i<=n;i++)
  {
    u+=x[i];  // x is a vector object which was previously declared
  }
}
\endexample
\noindent It is slightly less efficient to initialize x in this way. 

For the BORLAND\cplus\ compiler both libraries have been 
compiled to enable overlays.

\XX{overlays}{BORLAND \cplus\  compiler and}
\X{Borland}

\endchapter


\htmlnewfile 
%chapno=7
\mychapter{Input and Output}
\X{Input and Output}
\X{iostream}
\X{operators << >>}
\XX{vector}{I/O operations}
Input and output methods for the \AD\ class objects  
have been implemented by
overloading of the \cplus\ stream I/O operators {\tt <<}
and {\tt >>} for all \AD\ classes.
The so-called {\it iostream} package, introduced in 
AT\& T release 2.0 of \cplus\ is fully supported. 
For a detailed discussion of \cplus\ stream I/O you should
consult your \cplus\ manual (and the appropriate header files).
In this section, we
discuss extension of stream I/O to \AD\ class objects. 
Formatted input and output uses the standard \cplus\ stream classes, 
{\tt istream}, {\tt ostream}, {\tt ifstream}, and {\tt ofstream}. 
In addition, unformatted (``binary'') file I/O is implemented
through two additional stream classes {\tt uostream} and {\tt uistream}.

\section{Formatted Stream I/O}
\XX{Input and Output}{formatted}
\X{formatting} 
\XX{manipulators}{setw}
\XX{manipulators}{setfixex}
\XX{manipulators}{setscientific}
\XX{manipulators}{setprecision}
Variables of type \adtypes\ may be input or output using the
standard stream operators and formatted using the standard
manipulators. A two additional manipulators are provided to
simplify selection of either fixed point or ``scientific''
(e-format) notation.
The use of the manipulators {\tt setw()} (set width),
{\tt setprecision()}, {\tt setfixed}, {\tt setscientific}
and {\tt endl} is illustrated below
\beginexample
#include <iomanip.h>
#include <fvar.hpp>
main()
{
  int n = 30;
  int m = 5;
  gradient_structure gs;
  dvar_matrix x(1, n, 1, m);
  .
  .   // code to calculate x
  . 
  // e format with default width and precision
  cout << setscientific << x << endl;

  // fixed point format
  cout << setw(13) << setprecision(3) << setfixed << x << endl;
  .
  .
  .  
}
\endexample
\noindent 
The first output statement will produce thirty lines with five numbers
per line in ``scientific'' (e-format) separated by one or more blanks.
The second output statement will produce a table with 30 rows and 5 columns
with each member of {\tt x} printed with a fixed decimal place
with 3 significant figures in a field 13 bytes wide.
If the manipulators are omitted, the output will appear with the
default width, precision, and format specified by the compiler and with each
number separated by one spaces.

\goodbreak
The following standard manipulators are also supported:
\XX{manipulators}{setf}
\beginexample
setf(f, ios::floatfield);  // to select fixed or scientific
setf(b, ios::basefield);   // to base of number
setf(a, ios::adjustfield); // to select justification of fields
// where
//   f  may be ios::fixed or ios::scientific
//   b  may be ios::dec, ios::oct or ios::hex
//   a  may be ios::left, ios::right or ios::internal
\endexample

%Some compilers to not yet support {\tt iostream} package. The \AD\ system
%includes the overloaded function {\tt fform( )} which mimics the
%older style \cplus\ formatting function {\tt form(const char* ... )}.
%The function prototypes for {\tt fform( )} are:
%\beginexample
%char* fform(const char*, dvariable&);
%char* fform(const char*, dvar_vector&);
%char* fform(const char*, dvar_matrix&);
%char* fform(const char*, dvector&);
%char* fform(const char*, dmatrix&);
%\endexample
%\noindent The first argument is any valid C format string as described
%in the documentation for {\tt printf}. The second
%argument is the object to be printed.
%The function {\tt fform( )} is noticeably slower than the {\tt iostream} 
%package and the size of the object it can print is limited by
%by the size of an internal buffer. It offers no known advantages 
%over the {\tt iostream} package and is only provided to enable some
%control over the format of the output.
%{\bf There are incompatibilities between the new iostream package and the
%older stream package. Users are cautioned against trying to use both.}
 
\section{Error checking}
\X{I/O Error checking}
\X{operator !}
\XX{operator !}{I/O error checking and}
Although the stream operators {\tt >>} and {\tt <<} make it very easy to
perform input or output on the \AD\ container classes, it is important to
remember that they do not perform any error checking. While a complete 
discussion of error states for the stream classes is beyond
the scope of this manual a brief discussion of the use of the 
operator~{\tt !} for I/O error checking is included. If
{\tt infile} is a stream object then after any I/O operation
involving {\tt infile}, {\tt !infile} will be true if the operation
has failed and false if the operation has succeeded. 
For example:
\beginexample
dvector x(1,10);
ifstream infile("data"); // associate the file data with the ifstream 
                         // object infile
if (!infile) // Check if the file has been successfully opened
{
  cerr << "Error trying to open file data\n";  
  exit(1); // User has decided to stop on this error condition
}
                            
infile >> x ;            // read data into the dvector x;
if (!infile) // Check if the input operation was successful
{
  cerr << "Error trying to read dvector x from file data\n";  
  exit(1); // User has decided to stop on this error condition
}
\endexample
If you neglect to put in any error checking for I/O operations, your
program may well carry on blissfully even though the I/O operation has
failed. The symptoms of this failure can be quite confusing
and waste a lot of valuable program development time.

\section{Unformatted Stream I/O}
\XX{Input and Output}{unformatted}
\XX{Input and Output}{binary}
It is often useful to store intermediate calculations in an external
file. Unformatted (``binary'') file I/O is the most rapid and space 
efficient method for such operations and
is supported
through two additional stream classes {\tt uostream} and {\tt uistream}.
These two classes are derived classes from the
{\tt ofstream} and {\tt ifstream} classes in the {\it iostream} package. 

\section{An example of input and output for \AD\ classes}
\XX{Input and Output}{example}
The following listing illustrates the main features of the \AD\ system
I/O support. It also illustrates the use fill operations
and \cplus\ scoping rules to
automatically destroy some variables.
\bigskip
\listing{../../../examples/autodif/io_examp/io_examp.cpp}

\endchapter

\htmlnewfile 
%chapno=8
\mychapter{Temporary Files}
\X{gradfil1.tmp}
\X{gradfil2.tmp}
\X{cmpdiff.tmp}
\XX{temporary files}{{\tt gradfil1.tmp}}
\XX{temporary files}{{\tt gradfil2.tmp}}
\XX{temporary files}{{\tt cmpdiff.tmp}}
\XX{temporary files}{for saving derivative information}
The \AD\ system saves information required to calculate the
derivatives on a stack during the evaluation of the function. 
The amount of information depends on the number of variables for which
derivative information is required and on the specific nature of the
computations. It is difficult to predict generally how much information
will be generated. If the
length of this stack exceeds a preset amount, it is stored on the disk 
in one or two unformatted temporary gradient files, 
{\tt gradfil1.tmp} and {\tt gradfil2.tmp}. There is also
derivative information from ``precompiled'' derivative calculations
which is stored in the file {\tt cmpdiff.tmp}. 

\XX{RAM disk}{and temporary files}
\XX{temporary files}{use of the DOS environment string TMP}
\XX{temporary files}{use of the DOS environment string TMP1}
\XX{temporary files}{default directory for}
\XX{environment string}{specifying directory for temporary files}
The location of these files
is controlled by the DOS environment strings, {\tt TMP} and {\tt TMP1}.
The file {\tt gradfil1.tmp} will be created in the directory
indicated by {\tt TMP} while the file
{\tt cmpdiff.tmp} will be created in the directory indicated by 
 {\tt TMP1}. If either or both of these strings is undefined the 
current directory will be used to create the files.
Greatest speed will be achieved if the temporary files are located on a
RAM disk. Assuming that your RAM disk is drive {\tt E:}, type 
{\tt set tmp=e:} at the DOS prompt. If {\tt gradfil1.tmp}, becomes
too large for your RAM disk (or whatever device is specified by the
{\tt TMP} environment string), the rest of the file will be stored in
{\tt gradfil2.tmp} on 
the device referred to by the DOS environment string {\tt TMP1}. 
For example if you have typed {\tt set tmp=e:} and
{\tt set tmp1=c:\bs temp}, 
the first part of the gradient information will be stored on {\tt e:} 
(presumably a RAM disk) while the second part will
be stored in the directory {\ttt c:\bs temp}. 
If the directory does not exist, \AD\ will display
an error message and retire. If either the {\tt TMP} or {\tt TMP1} 
environment string does not
exist \AD\ will create {\tt gradfil1.tmp} or {\tt gradfil2.tmp} 
in the current directory.

\XX{temporary files}{effect of abnormal program termination on}
The temporary gradient files are normally deleted when objects of class 
{\ttt gradient_structure} go out of scope. Abnormal termination of
the program (such as by pressing {\tt Ctrl-Break} or use of the derivative
checker) will cause some temporary gradient files to remain
in your file system. These files may become quite large, several~Mb or more,
 so it is a good idea to remove them.
\endchapter

\htmlnewfile 
%chapno=9
\mychapter{Global Variables}
\section{Adjusting the \AD\ System Global Variables}

\XX{global variables}{access functions for}\X{gradient\_structure}
\X{controlling memory allocation}
The \AD\ system declares a number of global variables which
control memory allocation for the \AD\ structures.
These global variables can be changed by the use of
``access functions'' prior to declaring an object of type
{\tt gradient\_structure}.  
For most applications the
default values of these global variables should suffice. You will
not need to worry about them until you start building larger applications.
\beginexample
int NUM\_RETURN\_ARRAYS = 10;  
\endexample
\XX{NUM\_RETURN\_ARRAYS}{default value}  
\noindent This global variable determines the maximum allowable depth of 
nesting of functions which return \AD\ variable types (see the discussion about
the {\tt RETURN\_ARRAYS\_INCREMENT()} and 
{\tt RETURN\_ARRAYS\_DECREMENT()} instructions).
The default value of 10 will suffice for almost all applications and
most users will not need to concern themselves with this variable.
To increase the number to 15 you should put the instruction
\X{RETURN\_ARRAYS\_INCREMENT}\X{RETURN\_ARRAYS\_DECREMENT}
\XX{Error messages}{Overflow in RETURN_ARRAYS stack -- 
Increase NUM_RETURN_ARRAYS}
\beginexample
gradient\_structure::set\_NUM\_RETURN\_ARRAYS(15);  
\endexample
\noindent into your program before an object of type {\tt gradient\_structure}
is declared to manage the derivative calculations.
If you exceed the maximum allowabel depth of nesting the message
\beginexample
Overflow in RETURN_ARRAYS stack -- Increase NUM_RETURN_ARRAYS
\endexample
will appear. This message can also occur if you have put a
{\tt RETURN\_ARAYS\_INCREMENT} which is not matched by a
{\tt RETURN\_ARAYS\_DECREMENT} into one of your functions.

\beginexample
long int GRADSTACK\_BUFFER\_SIZE = 2200;  
\endexample
\X{GRADSTACK\_BUFFER\_SIZE}  
\X{CMPDIF\_BUFFER\_SIZE}  
\X{set\_CMPDIF\_BUFFER\_SIZE}  
\noindent 
This global variable determines the number of entries which are contained in
the buffer which contains the information necessary for caclulating 
derivatives. For historical reasons the actual amount of memory reserved for the buffer in bytes is equal to the value of
{\tt GRADSTACK\_BUFFER\_SIZE} 
multiplied by the size in bytes of an \AD\ structure,
{\tt grad\_stack\_entry}. For 16 bit DOS compilers the size of
{\tt grad\_stack\_entry} is 28 bytes. The default value of 
{\tt GRADSTACK\_BUFFER\_SIZE} is 2200 which reserves a buffer 
of~28*2200=61600~bytes.  Since for 16 bit DOS versions of \AD\ buffers 
are limited in size
to a maximum of 64K, 2200 is about the largest size which can be used.

If you make the value smaller (say 500) 
it will free up some memory (about 50K),
but the program will execute more slowly because it must store
data on the hard disk more often. This is a desperation move
in a situation where you must find some more memory for your program. 
To decrease the number to 500 you should put the instruction
\beginexample
gradient\_structure::set\_GRAD\_STACK\_BUFFER\_SIZE(500);  
\endexample
\noindent 
into your program before declaring a {\tt gradient\_structure} object.
\beginexample
long int CMPDIF_BUFFER_SIZE = 32000L;  
\endexample
This global variable determines the size in bytes of the buffer used to contain
the information generated by the ``precompiled'' derivative code.
To change this variable the instruction
\beginexample
static void gradient\_structure::set\_CMPDIF\_BUFFER\_SIZE(long int i);
\endexample
\noindent is used.

\bestbreak
\beginexample
int RETURN_ARRAYS_SIZE = 50;  
\endexample
\X{RETURN\_ARRAYS\_SIZE}  
\noindent This global variable controls the amount of complexity which one
line of arithmetic can have. The present default value should 
be large enough for any
conceivable purpose. It is recommended that the user leave it alone.

\beginexample
int MAX_NVAR_OFFSET = 200;  
\endexample
\X{MAX\_NVAR\_OFFSET}  
\XX{independent variables}{maximum number of}
\noindent This global variable determines the maximum number of independent
variables which can be used. 
In can be increased, for example, to 500 by including the instruction
\beginexample
\tt gradient\_structure::set\_MAX\_NVAR\_OFFSET(500);  
\endexample
\noindent into your program before declaring a {\tt gradient\_structure} object.
If you have more independent variables than the value of
{\tt MAX\_NVAR\_OFFSET} the error message 
\beginexample
You need to increase the global variable MAX_NVAR_OFFSET to xxx
\endexample
\noindent will appear and execution will be halted.
\XX{Error messages} {You need to increase the global variable 
 MAX\_NVAR\_OFFSET to xxx}
\beginexample
unsigned long ARRAY_MEMBLOCK_SIZE=100000L;  
\endexample
\X{ARRAY\_MEMBLOCK\_SIZE}  
\noindent
This  global variable determines the maximum amount of memory (in bytes)
available for \AD\ variable type container class objects ({\tt dvar\_vector,
dvar\_matrix $\ldots$ etc}). It represents the largest
single amount of memory which is required when an object of
type {\tt gradient\_structure} is declared. An \AD\ variable array object
requires about 8 bytes for each element of the array. A
{\tt dvar\_vector} object of size 500 will require slightly more than
4000 bytes of memory. A {\tt dvar\_matrix} object with 20~rows and 
30~columns would require 20 times 30 times 8 = 4800 bytes of memory.
To provide the maximum amount of memory for other purposes you can
reduce the size of 
{\tt ARRAY\_MEMBLOCK\_SIZE} until an error message
telling you to increase   {\tt ARRAY\_MEMBLOCK\_SIZE} appears. 
\XX{Error messages}{Need to increase ARRAY\_MEMBLOCK\_SIZE parameter}
\beginexample
unsigned MAX_DLINKS  = 1000; 
\endexample
\noindent The access function
\XX{dvariable}{controlling the maximum number of with set\_MAX\_DLINKS}
\XX{set\_MAX\_DLINKS}{setting the maximum number of dvariable objects}
\beginexample
gradient_structure::set_MAX_DLINKS(int i);
\endexample
\noindent determines the maximum number of {\tt dvariable} objects
which can be in scope at one time. The default value is 1000
which should be enough for most purposes. If the value is
exceeded the message 
\XX{Error messages}{ Need to increase the maximum number of dlinks}
\beginexample
Need to increase the maximum number of dlinks
\endexample
\noindent will appear and program execution will be terminated.

The complete list of access functions for the \AD\ system global variables is:
\XX{access functions}{for global objects}
\beginexample
static void set_NUM_RETURN_ARRAYS(int i);
static void set_ARRAY_MEMBLOCK_SIZE(unsigned long i);
static void set_GRADSTACK_BUFFER_SIZE(long int i);
static void set_MAX_NVAR_OFFSET(unsigned int i);
static void set_MAX_DLINKS(int i);
static void set_CMPDIF_BUFFER_SIZE(long int i);
\endexample
\noindent These functions may not be used while an object of
type {\tt gradient\_structure} is in scope. Doing so
will produce an error message and program execution will be halted.

\section{Stack Size}
\XX{stack}{allocating sufficient space for the}  
The size of the stack is usually controlled by a compiler-specific global
variable. The default stack size is generally too small for \AD.
It is recommended that the stack size be increased to 15000 or 20000 bytes.
The \AD\ library routines have been
compiled without the option of checking 
for stack overflow so that if you do not
reserve enough area for the stack, your program will probably crash
the system without any warning or simply behave in some strange manner. 
\XX{stack}{checking for overflow of the}

Stack size is set by inserting an appropriate line in the code prior
to the {\tt main()} block.
The following examples set the stack length to 15000 bytes.
Refer to your compiler documentation for more details.

\X{\_stklen} \X{Borland}
\XX{stack}{setting the stack size for Borland \string\cplus}  
For the Borland compiler, the appropriate syntax is:
\beginexample
extern unsigned _stklen = 15000; 
main()
{
  // ... your code
}
\endexample
\X{\_stack} \X{Zortech}
\XX{stack}{setting the stack size for Zortech \string\cplus}  
\noindent For the Zortech compiler, the appropriate syntax is:
\beginexample
unsigned int _stack = 15000; 
main()
{
  // ... your code
}
\endexample

\endchapter

\htmlnewfile 
%chapno=10
\mychapter{The \AD\ function minimizing routines}
\XX{minimization of functions}{quasi--Newton method}
\XX{minimization of functions}{conjugate gradient method}
\X{conjugate gradient method}
\X{quasi--Newton method}
\AD\ is supplied with both a quasi-Newton function minimizer
and a conjugate gradient function minimizer.
A quasi-Newton function minimizer
is an effective method for minimizing smooth functions of up to one 
or two hundred parameters. At some point, the calculations involved in
updating the inverse Hessian approximation (on the order
$n^3$) and the memory requirements for
the inverse Hessian approximation ($n(n+1)/2$ floating point elements), 
where $n$ is the number of parameters, become prohibitive
and the conjugate gradient method, supplied with \AD, should
be used.  

\XX{minimization of functions}{user interface}
We have provided the user
with two types of interfaces to the function minimizer routines. The first approach provides the user with the maximum amount of encapsulation of the
minimization routine. It effectively allows the user to say ``Here is my
function. Minimize it.'' The disadvantage of this approach is some loss of
flexibility in managing memory and controlling the minimization
process. To provide maximum flexibility we also provide a lower level interface
which allows the user to exercise complete control over the minimization 
process.

\X{minimize function} 
The fully encapsulated minimization procedure is invoked by calling the
function {\tt minimize}. 
The following code invokes the {\tt minimize} routine
to find the minimum for the user's function {\tt userfun}.
\bigskip
\listing{../../../examples/autodif/minimize/minimize.cpp}

\XX{independent\_variables}{declaration of}
\XX{independent\_variables}{initialization of}
In this example, it is the user's responsibility to determine the
number of independent variables {\tt nvar} and to declare the array
of type {\tt independent\_variables} which holds the independent variables.
 The elements of the
{\tt independent\_variables} array {\tt x} are initialized to 0.  If
other initial values are desired, the user must initialize the
array as well. An example of a typical simple function which could
be minimized by this code is the code for the function
$$\big(x[1]-1\big)^2 + \sum_{i=1}^{i<n} \big(x[i+1]-x[i]\big)^2$$
whose minimum is at the point $(1,1,\ldots,1)$.
\bigbreak
\listing{../../../examples/autodif/minimize/userfun.cpp}

\X{dvar\_vector} \X{dvariable} 
While the vector of independent variables is declared
to be an object of type {\tt independent\_variables} in the main 
routine it must be declared to be an object of type {\tt dvar\_vector}
in {\tt userfun}. This is the \AD\  ``array'' type for which
derivative 
information is calculated. It corresponds to an array or vector
of doubles. Also, the quantities {\tt z} and {\tt tmp} which correspond to
objects of type {\tt double} for which derivative information is required
must be declared to be of type {\tt dvariable}. The function {\tt value}
takes a {\tt dvariable} as its argument and returns a {\tt double}
which is the value of the {\tt dvariable}.

This simple implementation of the function minimizer 
is useful for the beginning user who 
doesn't wish to get too involved in the intricacies of \AD\ or
who wants a quick solution to a minimization problem. What are its
disadvantages? 
The main disadvantage is that the type of arguments taken by the function
{\tt userfun} are limited to a predefined set. In the present case the
function {\tt userfun} is declared to be a function which takes
one argument of type {\tt dvar\_vector\&}. The line of code 
\beginexample
double minimum_value= minimize(x,userfun);
\endexample
\noindent causes a pointer to the function {\tt userfun} to be passed
to the function {\tt minimize} (when a function is named without the
accompanying parentheses a pointer to the function is generated).
Thus the second argument of {\tt minimize} is a pointer to
a function which takes one argument of type {\tt dvar\_vector\&}
and returns a {\tt double}.


To enable the user to retain maximum control over memory management
and the minimization process,
the function minimizer and the 
function being minimized should be put on the same level in the program.
 This is accomplished
by putting them into a loop. The general form of this control
structure is  
\beginexample
double f;
independent_variables x(1,nvar);
dvector g(1,nvar);
fmm fmc(nvar);  // Use this declaration for the quasi-Newton function minimizer
// fmmc fmc(nvar);  // Use this declaration for the conjugate-gradient function 
                    //minimizer
gradient_structure gs;  
while (fmc.ireturn >=0)
{
  fmc.fmin(f,x,g);
  if (fmc.ireturn > 0)
  { 
    f=userfun(x,...) // The users function
    gradcalc(nvar,g);
  }
}
\endexample
The {\tt dvector} objects {\tt g} and {\tt x} must have
minimum valid index equal to 1 and maximum valid index equal to the
number of independent variables in the minimization.

\X{BEGIN\_MINIMIZATION}
\X{END\_MINIMIZATION}
Since many parts of the code sequence are independent of the 
particular application the code segments for calling the function minimizer and
evaluating the derivatives has been made into a macro. Using the
macros, {\tt BEGIN\_MINIMIZATION(nvar,f,x,g,fmc)} and 
{\tt END\_MINIMIZATION(nvar,g)} 
the above code can be written as
\beginexample
double f;
independent_variables x(1,nvar);
dvector g(1,nvar);
fmm fmc(nvar); // Use this declaration for the quasi-Newton function minimizer
// fmmc fmc(nvar);  // Use this declaration for the conjugate-gradient function 
                    //minimizer
BEGIN_MINIMIZATION(nvar,f,x,g,fmc) 
  f=userfun(x,...) // The users function
END_MINIMIZATION(nvar,g) 
\endexample

The macro {\tt BEGIN\_MINIMIZATION} takes five arguments, the number of
independent variables, the name of the object of type {\tt double}
which will contain the value of the function being minimized,
the independent variables, the gradient vector, and the name of 
the function minimizing control structure.
The macro, {\tt END\_MINIMIZATION} takes two arguments, the number of
independent variables and the gradient vector. 
\XX{BEGIN\_MINIMIZATION}{number of arguments} 
\XX{END\_MINIMIZATION}{number of arguments} 

\XX{fmm class}{controlling minimization with}
\XX{fmmc class}{controlling minimization with}
Notice that even when using the macros, it is the responsibility of the
user to declare 
the name of the object of type {\tt double}, {\tt f},
which will contain the value of the function being minimized, the 
{\tt independent\_variables x}, and the {\tt dvector g} as well
as the {\tt fmm} object {\tt fmc}.
The quasi-Newton function minimizer {\tt fmin} is 
controlled by and is a member function 
of the class {\tt fmm}. The constructor for the
class {\tt fmm} takes {\tt nvar}, the number of variables,
as an argument and creates the structure, {\tt fmc} which controls 
the minimization.  The user can control the operation of the 
function minimizer by changing certain values of the {\tt fmm}
structure {\tt fmc}.
The conjugate gradient function minimizer {\tt fmin} is
a member function of the class {\tt fmmc}.
The only difference in calling the quasi-Newton or conjugate-gradient 
function minimizers is in the declaration of an
object of type {\tt fmm} for the quasi-Newton method and
of type {\tt fmmc} for the conjugate-gradient method.
The public members of this structure are given below. 
\XX{fmm class}{maxfn}
\XX{fmm class}{min\_improve}
\XX{fmm class}{crit}
\XX{fmm class}{imax}
\XX{fmm class}{iprint}
\XX{fmm class}{scroll\_flag}
\XX{fmm class}{ireturn}

\beginexample
class fmm
{
public:
                        // control variables:

  long   maxfn;        // maximum number of function evaluations
                       // default value: maxfn = 500
  double min_improve;  // if the function doesn't decrease by at least this
                       // amount over 10 iterations then quit.
                       // default value: min_improve=1.e-6
  double crit;         // gradient convergence criterion; fmin terminates if
                       //   maximum absolute value of the
                       //   gradient component is less than crit.
                       //   default value: crit = 1e-4
  long   imax;         // maximum number of function evaluations in linear
                       // search without improving the current estimate
                       // before fmin terminates  default value : imax = 30
  long   iprint;       // number of iterations between screen displays of
                       // current parameter estimates default value: iprint = 1
  long   scroll_flag;  // set to 1 to scroll display on screen; 
                       //        0 to overwrite screen 
                       // default value: scroll_flag = 1
  int    ireturn;      // used to control the loop containing fmin and the
                       // user's function 
} 
\endexample
There are five criteria by which the function minimizer may decide that
convergence has been achieved.
\X{convergence criteria for function minimizer}
\XX{function minimizer}{convergence criteria for} 
{\parindent=30truept
\begin{description}
\myitem{1.} If the magnitude of all the gradient components is less
than {\tt crit} (default value 1.e-4).

\myitem{2.} If the maximum number of function evaluations, {\tt maxfn} is exceeded
(default value 500)

\myitem{3.} If the function value does not decrease by more than
{\tt min\_improve} over 10 iterations (default value 1.e-6).

\myitem{4.} If the function minimizer tries to evaluate the function
more than {\tt imax} times within one linear search. (default value 30)

\myitem{5.} If the function minimizer detects that the user has pressed 
the Q or q key. 
\end{description}
\par }
Here is an example of the use of these switches.
\beginexample
dvector g(1,nvar);
dvector x(1,nvar);
double f;
fmm fmc(nvar);
fmc.maxfn=1000;       // Maximum number of function evaluations is 1000
fmc.iprint=10;        // Current best estimate will be printed out every
                      // ten iterations
fmc.crit=.01;         // Convergence if derivatives are all smaller in 
                      // magnitude than .01
fmc.min_improve=0.0;  // This convergence criterion won't be used
fmc.imax=0;           // This convergence criterion won't be used
gradient_structure gs;  
while (fmc.ireturn >=0)
{
  fmc.fmin(f,x,g);
  if (fmc.ireturn > 0)
  { 
    f=userfun(x,...) // The users function
    gradcalc(nvar,g);
  }
}
\endexample



The function prototype for {\tt fmin} is:
\beginexample
void fmm::fmin(double& f, independent_variables & x, dvector & g);
/*---------------------------------------------------------------------------*
 * Purpose: quasi-Newton numerical minimization                              *
 *  f       : current value of function to be minimized                      *
 *  x       : current values of the variables over which the minimization    *
 *            occurs; instance of class independent variables with elements  *
 *            from 1 to n                                                    *
 *  g       : current value of the gradient as returned by the function      *
 *            gradcalc( ... ); instance of class dvector with elements from  *
 *            1 to n                                                         *
 *---------------------------------------------------------------------------*/ 
\endexample


The execution of {\tt fmin} may be terminated by pressing the `{\tt Q}'
key on the keyboard. There may be a slight delay between the keypress
and the termination of {\tt fmin}.
\section{Putting bounds on parameters}

It is often necessary to put bounds on the values that a parameter
can take. For example if the expression $1/y$ is to be evaluated
in a program an arithmetic error will occur if $y=0$.
Similarly $e^x$ will probably produce an error if $x=100000$.

\X{boundp}
\XX{bounding functions}{boundp} 
The function {\tt boundp} is supplied with \AD\ to simplify the
placing of bounds on parameters.

\beginexample
dvariable boundp( dvariable x, double fmin, double fmax, dvariable& fpen)
{
  dvariable t,y;
  t=fmin + (fmax-fmin)*(sin(x*1.570795)+1)/2;
  if (x < 0.00001)
  {
    fpen+=(x-0.00001)*(x-0.00001);
  }
  if (x > 0.9999)
  {
    fpen+=(x-0.9999)*(x-0.9999);
  }
  if (x < -1)
  {
    fpen+=1000*(x+1)*(x+1);
  }
  if (x > 1)
  {
    fpen+=1000*(x-1)*(x-1);
  }
  return(t);
}
\endexample 
The function {\tt boundp} constrains the variable to lie
between the values {\tt fmin} and {\tt fmax}. Both the values
{\tt fmin} and {\tt fmax} can actually be attained so if you
want a variable to be greater than 0 the code.
\beginexample
y=boundp(x,0.0,1.0,fpen); 
\endexample
\noindent 
will not do the job. Instead you can accomplish this in one of two ways.
You can use a quadratic transformation if you really don't
care how large x can become.
\beginexample
y=x*x; 
\endexample
\noindent You can make a decision about what is meant by ``close to~0''
for your application, in this case 0.00001.
\beginexample
y=boundp(x,0.00001,1.0,fpen); 
\endexample

\XX{penalty function}{use wish bounding functions}
The penalty function {\tt fpen} is included in the {\tt boundp}
routine because otherwise the function would have a zero 
derivative at its lower and upper bounds. This causes the
derivative with respect to the bounded parameter to 
be equal to 0. The result is that if the minimization routine
is started with a parameter value already at the upper or
lower bound then this parameter will never be changed.
The penalty shifts the position of the  critical point
slightly within the bounds. Of course if the initial parameter
estimate happens to coincide exactly with this critical point the
derivative will be 0, but this is very unlikely.

If you neglect to add the penalty function to the function
being minimized the minimization routine will still work.
The only problem is that you run the risk of one of the
bounded parameters getting ``stuck'' at its upper or lower bound.

\XX{penalty function}{getting the sign correct}
The penalty should be added to the function which you wish
to minimize. If you wish to maximize the function (in which case
you are going to change the sign of the function before passing it to
the function minimizer) then you should subtract the penalty.
{\bf If you get the sign of the penalty wrong the routines will not
work properly.}

\X{boundpin}
\XX{bounding functions}{boundpin}
\XX{bounding functions}{inverse function for}
To begin the minimization you generally know the initial value of the
constrained parameter~{\tt y}, and you need to find the value of~{\tt x}
which will produce this value so that you can initialize the 
proper component of the vector of {\tt independent\_variables}
with it.  The function {\tt boundpin} which is the
inverse function of {\tt boundp} will yield the
desired value for~{\tt x}.  
\beginexample
// double boundpin(double y, double fmin, double fmax)
x=boundpin(y,fmin,fmax);
\endexample

\section{The derivative checker}
\XX{Derivative checker}{checking correctness of derivatives with}
It may happen that the derivatives being passed to fmin become incorrect
for some values of the parameters. Two possible reasons for this are,
an error in the internal derivative calculations of \AD\ itself or,
the use of a user's function which is not differentiable at the point  
in question. The former cause becomes increasingly less likely as
experience with \AD\ continues and any remaining bugs are removed. The latter
can occur because when constructing a large complicated function it is always
possible (and not uncommon in our experience) to add some non-differentiable 
components as an oversight.

If the function becomes nondifferentiable at some point or the
derivative calculations are incorrect, the function minimizer will
begin to experience difficulty in making progress.
To help identify derivative errors or cases of non differentiability 
when they occur
The \AD\ system includes a interactive derivative checker which compares the
value of the derivatives calculated by {\tt gradcalc()} with 
finite difference approximations. The derivative checker 
can be invoked while the
function minimizer is in progress by pressing the `{\tt C}' key on the keyboard.
\XX{Derivative checker}{invoking from the function minimizer}

The derivative checker operates in two modes. In 
the default mode, the derivative is checked for one independent variable 
at a time. In this mode you will be
asked to supply the index of the variable for which you want to check
the derivative. 
Alternatively, if you want to check all derivatives, enter $0$ when
prompted for the index of the variable to check. Next
you will be asked to enter the step size to use for the finite difference
approximation. 
Enter a number which will alter the value of the of the variables
to be checked in the fifth or sixth significant digit. 
The derivative checker will
then tabulate the values of the function, the variable being checked,
the analytical derivative (from {\tt gradcalc()}), the finite difference
approximation, and the index of the variable being checked.
The derivative checker may be interrupted by pressing the `{\tt X}'
key on the keyboard. On completion, the derivative checker exits to
the operating system. 

For many models the finite approximations to the derivatives may
change radically as the step size used in the derivative approximation 
is changed. How then can one be sure that the analytical derivatives 
are correct or incorrect. The following ``rule of thumb'' seems to
work fairly well. First identify a gradient component which seems
to be incorrect. Pick this component as the independent-variable
whose derivative is to be checked. Evaluate the finite difference
approximation for different step sizes, with the step sizes
differing by a factor of about 3. Typical step sizes might be
1.e-4, 3.e-5, 1.e-5, 3.e-6, $\ldots$. Now if the finite difference
approximations stay almost the same for two consecutive step sizes,
but are very different from the analytical approximations, this is
a good indication that the analytical derivatives are incorrect.

 It is wise to use the derivative checker
frequently during the development of an application and after any change
the the calculation of the objective function.

\XX{Derivative checker}{invoking explicitly}
The derivative checker may also be invoked explicitly. Its function
prototype is:
\beginexample
void derch(double & f, independent_variables & x, dvector & g,
           const int n, int & ireturn);
/*---------------------------------------------------------------------------*
 * Purpose: compare analytical derivative calculation with finite difference *
 *          approximation                                                    *
 * Input:                                                                    *
 *  f       : current value of function                                      *
 *  x       : current values of the parameters of the function; instance of  *
 *            class independent_variables with elements ranging from 1 to n  *
 *  g       : current value of the gradient as returned by the function      *
 *            gradcalc( ... ); instance of class dvector with elements from  *
 *            1 to n                                                         *
 *  n       : number of variables                                            *
 *  ireturn : return status; set ireturn = 3 prior to first call to derch()  *
 * Output:                                                                   *
 *  ireturn : return status; ireturn = 4 or 5 on return from derch()         *
 *---------------------------------------------------------------------------*/           
\endexample

The following code fragment will display 
the values of the function, the variable being checked,
the analytical derivative (from {\tt gradcalc()}), the finite difference
approximation, and the index of the variable being checked.

\beginexample
#include <math.h>
#include <fvar.hpp>

// prototype for user supplied code to compute f as a function of x
void f_comp( double & f, independent_variables & x, ... ); 

void main()
{
  int n = ...;        // number of variables in function    
  independent_variables x(n);
  set_x( ... );      //  user supplied code to set the values of x
  dvector g(1,n);
  gradient_structure gs;  
  double f = 0;
  int derch_return = 3; 

  f_comp(f, x, ... );   // first call to get value of function for current x
  while (derch_return >=3)
  { 
    derch(f, x, g, n, derch_return);
    {
      f_comp(f, x, ... );    // subsequent calls for with values of x
                             // varied by the step size within derch()
      gradcalc(n, g);
    }
  }
}
\endexample
\endchapter

\htmlnewfile 
%chapno=11
\mychapter{Statistical functions}

The principal purpose of \AD\ is for constructing computer programs
that use mathematical derivatives. One such application
is in nonlinear  parameter estimation, that is the construction
of nonlinear statistical models. 
It is convenient to use simulations which produce 
data with known statistical properties in order
to analyze the behaviour of such models.
To facilitate this
process \AD\ is furnished with a small number of statistical
functions which can be applied to both constant and variable vector objects.
\X{statistical functions}
\X{simulation models}
\X{nonlinear parameter estimation}

\section{Filling vectors with random numbers}

\XX{random numbers}{filling vectors with}
\XX{filling vectors} {with random numbers}
In this section, a uniformly distributed random number is assumed to have
a uniform distribution on $[0,1]$. A normally distributed random number
is assumed to have mean $0$ and variance $1$.
A binomially distributed random number is assumed to have a parameter $p$
where $1$ is returned with probability $p$ and $0$ is returned with
probability $1-p$.
The {\tt long~int\&~n} is the random number seed specified by the user.

\X{fill_randu}
\X{fill_randn}
\X{fill_randbi}
\X{fill_multinomial}
\beginexample
dvector::fill_randu(long int& n) // Fill a dvector from a uniform 
                                 // distribution on 0,1 
dvector::fill_randn(long int& n) // Fill a dvector from a standard normal 
                                 // distribution 
dvector::fill_randbi(int& n,double p) // Fill a dvector from a binomial 
                        //distribution A 1 is returned with probability p
                        // A 0 is returned with probability 1-p
dvector::fill_multinomial(int& n,dvector p)// Fill a dvector from a multinomial 
             //distribution The distribution is determined by the dvector p.
            // An i is returned with probability p(i) if the components of p
            // sum to 1. Otherwise p is normalized so that its components
            // sum to 1.  
ivector::fill_multinomial(int& n,dvector p)// Fill a dvector from a multinomial 
             //distribution The distribution is determined by the dvector p.
            // An i is returned with probability p(i) if the components of p
            // sum to 1. Otherwise p is normalized so that its components
            // sum to 1.  
\endexample
The {\tt fill\_multinomial} routine expects to receive a {\tt dvector}
of parameters $P=(p_1,\ldots,p_n)$. The routine fills the the elements of
a {\tt dvector} with integers whose values range from $1$ to $n$ 
where $i$ occurs with probability $p_i$.
If the components of $P$ do not sum to $1$ they will be normalized
so that the components do sum to $1$.

As an example of how to use these functions, consider the following
code segment which generates random numbers from a normal distribution
which is contaminated with a small proportion of large errors.
This kind of distribution is useful for studying ``robust'' parameter 
estimation procedures. Robust estimation procedures are intended to
be insensitive to such deviations from the main distribution.
\XX {random number generator}{example of the use of}\X{robust estimation} 

Assume that the main distribution has mean~20 and standard deviation~2,
while the contaminating distribution has mean~20 and standard 
deviation~6. We want the probability that the random numbers belong to
the main distribution to be 0.90 and the probability that they 
belongs to the contaminating distribution to be 0.10. We shall
generate a random vector of 100 such random numbers. 
\beginexample
// ...
dvector rn(1,100);
{
  dvector w(1,100);
  rn.fill_randn(212); // fills rn with numbers drawn from a normal distribution 
  w.fill_randbi(1521,0.90); // fills w with numbers drawn from a binomial
                            // distribution   
  for (int i=1; i<=100;i++)
  {
    if (w(i)==1)  // This condition will occur with probability 0.90
    {
      rn(i)=2*rn(i)+20;  // standard deviation 2 and mean 20
    }
    else// This condition will occur with probability 0.10
    {
      rn(i)=6*rn(i)+20;  // standard deviation 6 and mean 20
    }
  }
} // Now the dvector w goes out of scope and the memory it used is released
// ...
\endexample

\section{Generating a sample from a mixture of n normal distributions}
\X{mixture of normal distributions}  
\XX{multinomial fill}{of an ivector object}
This example extends the previous example to a mixture of more than two
distributions. Two interesting aspects of this example are the use of
the multinomial fill of an {\tt ivector} object and the use of the resulting
{\tt ivector} object in a general selection of the form {\tt dvector(ivector)}.
\X{operator ()}
\XX{vector}{extracting a subvector}
\X{extracting a subvector}
\XX{vector}{function call () to extract subvector}

\beginexample
int nsample;  // nsample is the sample size to be generated
int n;        // n is the number of groups in the mixture
dvector p(1,ngroups); // p contains the mixture proportions;
dvector mean(1,ngroups); // mean contains the mean for each group 
dvector sd(1,ngroups); // sd contains the standard deviations for each group 
dvector sample(1,nsample); // sample will contain the simulated data
ivector choose(1,nsample); // this will determine which group in the mixture
                           // each observation came from
//...
// somehow get the mixture proportions into p, the means into mean,
// and the standard deviations into sd
// ...
choose.fill_multinomial(1011,p); 

sample.fill_randn(211); // fill the sample with standard normal deviates
sample=elem_prod(sd(choose),sample)+mean(choose); //This
            // creates the mixture of normal distributions
\endexample
Notice that the $i$'th element of {\tt choose} will be equal to $j$ with
probability $p_j$ so that each element in sample will have
a probability $p_j$ of having mean(j) as its mean and sd(j) as its standard
deviation. This means that each element of {\tt sample} will
have a probability $p_j$ of belonging to the $j$'th group of normal
distributions. It follows that the elements of {\tt sample} have
the desired mixture distribution.
\section{Functions which provide summary statistics}
\XX{mean}{mean of a vector object}
\XX{std\_dev}{standard deviation a vector object}
The functions
\beginexample
number mean(vector_object);
number std_dev(vector_object);
\endexample
\noindent return the mean and standard deviation of a {\tt dvector} or 
{\tt dvar\_vector}.
\endchapter


\htmlnewfile 
%chapno=12
\def\mynarrower{\advance\leftskip by 20truept\advance\rightskip by 20truept}
\def\ls{least squares }
\def\ei{\epsilon_i}
\def\mm{least median of squared residuals }
\mychapter{Robust Nonlinear Regression}
Many authors have indicated that the usual method of least 
squares estimation of parameters is too sensitive to the existence of a 
small number of large errors, so-called outliers, in the data being analyzed.
Although an extensive theory of robust estimators  which are not so 
sensitive to 
such outliers has developed, many people still use standard statistical
packages based on least squares estimation for statistical problems
such as multiple regression. With the \AD\ system it is a simple matter
to use superior
robust regression estimation procedures on both linear and
nonlinear statistical models.
We should emphasize that we are not experts in the field of robust regression
and we make no claims of optimality or near optimality for these methods.
 The methods in this chapter seem to work well and they
have been designed so that they can be easily extended to large 
nonlinear regression problems with thousands of observations
and at least hundreds and perhaps thousands of parameters. 
(In 16 bit DOS versions of \AD\ the 64K barrier may limit the
size of the problem which can be considered).
\XX{least squares}{parameter estimation}\XX{robust estimation}{in nonlinear parameter estimation}
\XX{parameter estimation}{least squares}
\XX{parameter estimation}{robust nonlinear}

Our purpose in this chapter is to present an estimation procedure
which performs almost as well as the usual least squares estimator
when the model errors are normally distributed and which performs much
better than least squares when the model errors contain some large
outliers. With such an estimator you don't have to worry
about whether the use of robust techniques is justified. If
your model errors are exactly normally distributed you have lost
very little and if they are not normally distributed you may 
have gained a lot.

The reader who is not interested in
the derivation of the formulas in this 
section can skip to section \number\mychapno.2 where the methods for
using the routines are discussed. 
The necessary functions are supplied with the \AD\ libraries. 

\section{Statistical Theory for robust regression}

Assume that we have $n$ data points, $Y_i$, $i=1,\ldots,n$. 
The $Y_i$ are assumed to have been produced by the process
$$\myeq{Y_i=f({\mbox{\boldmath $x$}}_i;\Theta)+\epsilon_i}\eqno(\number\mychapno.0)$$
where for each $i$, ${\mbox{\boldmath $x$}}_i$ is a known vector, $\Theta$
is an unknown vector of parameters, and the~$\epsilon_i$ are
random variables. We want to get estimates for the
parameters~$\Theta$. Some common special cases are:
\bestbreak
{\mynarrower\openup -1truept
Multiple regression:
\smallskip
for each $i$,
${\mbox{\boldmath $x$}}_i =(x_{i1},\ldots,x_{ir})$ is an $r$-vector and
$\Theta =(\theta_{0},\theta_{1},\ldots,\theta_{r})$ is an $r+1$-vector. The
functional relationship is
$$f({\mbox{\boldmath $x$}}_i;\Theta)=\theta_0+x_{i1}\theta_1+x_{i2}\theta_2+\ldots+x_{ir}\theta_r$$
\smallskip
Polynomial regression:
\smallskip
for each $i$, ${\mbox{\boldmath $x$}}_i =(x_i)$  is a $1$-vector and
$\Theta =(\theta_{0},\theta_{1},\ldots,\theta_{r})$ is an $r+1$-vector. The
functional relationship is
$$f({\mbox{\boldmath $x$}}_i;\Theta)=\theta_0+x_i\theta_1+x_i^2\theta_2+\ldots+x_i^r\theta_r$$
\smallskip
A nonlinear model:
\smallskip
For  a simple nonlinear model we shall condsider the case where ${\mbox{\boldmath $x$}}_i$
is a $1$-dimensional vector and $\Theta$ is a $1$-dimensional vector
whose single component is denoted by $b$. The functional
relationship is 
$$\myeq{f({\mbox{\boldmath $x$}}_i;b)=bx_i+b^2}\eqno(\number\mychapno.1)$$
\smallskip
}

The form of the estimation procedure which should be
used to estimate the parameter 
vector~$\Theta$  depends on the assumptions which are made about
the nature of the random variables~$\epsilon_i$. Assume that the~$\epsilon_i$ 
are independent identically distributed random variables
with a probability density function~$\phi(u;\Lambda)$ where~$\Lambda$
is a vector of parameters whose components may be known or unknown. 

The joint probability density function for the $(Y_1,\ldots,Y_n)$ is given by
$$\myeq{\prod_{i=1}^n \phi\big(Y_i-f({\mbox{\boldmath $x$}}_i;\Theta);\Lambda\big)}
\eqno(\number\mychapno.2)$$
Taking the logarithm of $(\number\mychapno.2)$ gives the expression
$$\myeq{\sum_{i=1}^n \log\Big(\phi\big(Y_i-f({\mbox{\boldmath $x$}}_i;\Theta);\Lambda\big)\,\Big)}
\eqno(\number\mychapno.3)$$
Considered as a function of the parameters $\Theta$, and $\Lambda$ expression $(\number\mychapno.3)$
is known as the log-likelihood function. The maximum likelihood estimates 
for the parameters $\Theta$ and $\Lambda$ is the value of parameters which 
maximizes the likelihood, or what is equivalent, the log-likelihood function.

If we assume that the $\epsilon_i$ are normally distributed with
mean~$0$ and standard deviation~$\sigma$ then
$$\phi(u;\Lambda)= {1\over\sqrt{2\pi}\sigma} 
    \exp\big(-{u^2\over 2\sigma^2}\big)$$
\noindent and the log-likelihood function becomes
$$\myeq{-n\log(\sigma)-{1\over2\sigma^2}\sum_{i=1}^n\big(Y_i-f({\mbox{\boldmath $x$}}_i;\Theta)\big)^2}
             \eqno(\number\mychapno.4)$$
In this case $\Lambda=(\sigma)$ is a 1 dimensional vector.

The maximum likelihood estimates estimates
$\widehat\Theta$ for the parameters $\Theta$ which maximize 
$(\number\mychapno.4)$ are the same as the those values of the parameters
which minimize
$$\myeq{\sum_{i=1}^n\big(Y_i-f({\mbox{\boldmath $x$}}_i;\Theta)\big)^2} \eqno(\number\mychapno.5)$$
Since these parameter estimates minimize the 
squared difference between the predicted and
observed values they are known as the least-squares estimates. 
For normally distributed model errors the maximum likelihood estimates
and the least squares estimates coincide. Since the expression
$(\number\mychapno.5)$ does not depend on~$\sigma$, The maximum-likelihood
estimate for~$\Theta$ does not depend on the estimate for~$\sigma$.
This is a special property shared by the normal distribution and
some other distributions. It is does not hold in general.
The maximum likelihood 
estimate~$\hat\sigma$ for the standard deviation~$\sigma$ is given by
$$\myeq{\hat\sigma=\sqrt{{1\over n}\sum_{i=1}^n\big(Y_i-f({\mbox{\boldmath $x$}}_i;\Theta)\big)^2}}
  \eqno(\number\mychapno.6)$$
\XX{maximum likelihood estimation}{for robust regression}

The reason that the normal least squares estimation is so sensitive to
outliers is that the normal distribution has a small ``tail'',
so that the probability of occurrence of an event drops
off very quickly as one moves away from the mean a distance of a few
standard deviations.  One approach to develop more robust
estimation methods is to use distributions with fatter tails.
A distribution with a fatter tail than the normal distribution
is the double exponential distribution whose probability
density function (assuming the mean is 0) is given by
\X{two-sided exponential distribution}
$$\myeq{{1\over2\sigma} 
    \exp\big|-{u\over\sigma}\big|}\eqno(\number\mychapno.7)$$
The corresponding log-likelihood function is
$$\myeq{-n\log(\sigma)-{1\over\sigma}\sum_{i=1}^n\big|Y_i-f({\mbox{\boldmath $x$}}_i;\Theta)\big|}
             \eqno(\number\mychapno.8)$$
The values of the parameters $\Theta$ which maximize
$(\number\mychapno.8)$ are the same as the values which minimize the sum
$$\myeq{\sum_{i=1}^n\big|Y_i-f({\mbox{\boldmath $x$}}_i;\Theta)\big|}\eqno(\number\mychapno.9)$$
Since $(\number\mychapno.9)$ is the total absolute deviation between
the predicted and observed values these estimates are
known as the minimum total absolute deviation estimates for the
parameters.

While the minimum total absolute deviation estimates  are fairly
insensitive to a small number of large  model errors, they
do not perform very well when the model errors are normally distributed.
To develop a maximum likelihood estimation
procedure which is less sensitive to outliers, but which
also performs well when the model errors are normally distributed,
 we have employed a distribution which is a 
mixture of a normal
distribution and another distribution with  a fatter tail.
The probability density functions of the fat-tailed distribution is given by 
$$\myeq{{\sqrt{2}\over\pi\sigma e}\Big\{1+{x^4\over{(\sigma e)}^4}\Big\}^{-1}}
     \eqno(\number\mychapno.10)$$
The extra parameter $e$ has been added to the log-likelihood function
$(\number\mychapno.10)$ to adjust the ``spread'' of the
fat-tailed distribution. The value of $e$ has been set equal to~3.
A mixture of two random variables with mixing proportion $p$ is
formed by choosing one of the random variables with probability~$p$
and the other random variable with probability~$1-p$. If the
probability density functions of the components of the mixture
are $\psi(u)$ and $\chi(u)$ then the probability density function
of the mixture is $p\psi(u)+(1-p)\chi(u)$.
We shall assume that the model errors $\epsilon_i$ are a mixture 
of the distribution $(\number\mychapno.10)$ with probability~0.05
and a normal distribution with probability~0.95.
The quantity $0.05$ can be interpreted as the
proportion of observations contaminated by large errors.  

The corresponding log-likelihood function for the observations is 
$$\myeq{-n\ln(\sigma)+\sum_{i=1}^n \log\bigg[
    (1-p)\exp\bigg\{-{\big(Y_i-f({\mbox{\boldmath $x$}}_i;\Theta)\big)^2\over 2\sigma^2}\bigg\}+
    {2p\over \sqrt{\pi} e}\bigg\{1+{\big(Y_i-f({\mbox{\boldmath $x$}}_i;\Theta)\big)^4
     \over (e\sigma)^4}\bigg\}^{-1}\bigg]}
\eqno(\number\mychapno.11)$$

The parameter estimates $\widetilde\Theta$ and $\tilde\sigma$ for $\Theta$
and $\sigma$ obtained by maximizing $(\number\mychapno.5)$ will
be referred to as the robust-mixture estimator for the parameters. 
For the robust-mixture estimator the estimate for $\Theta$
is not independent of the estimate
for $\sigma$, so that to estimate the parameters the entire
expression $(\number\mychapno.11)$ must be maximized for $\Theta$ and $\sigma$ simultaneously.

While it is true that the robust parameter estimates are the values
of the parameters which maximize $(\number\mychapno.11)$ it is possible to
reformulate the problem slightly  so that it is posed in terms of
parameters which are more easily estimated. Such reparameterizations are often
useful in nonlinear optimization. 
Consider the log-likelihood function $(\number\mychapno.11)$. If~$p$
is set equal to~0 the log-likelihood reduces to the normal
case, expression, $(\number\mychapno.4)$ 
so that the value of the parameter~$\sigma$ which would
maximize $(\number\mychapno.11)$ 
when~$p=0$ is~$\hat\sigma$ given by $(\number\mychapno.6)$.

Define a new parameter $\alpha$ by $\alpha=\sigma/\hat\sigma$ so that
  $\sigma=\alpha\hat\sigma$. 
Replace the parameter~$\sigma$ in $(\number\mychapno.11)$ by
the parameter~$\alpha\hat\sigma$. The log-likelihood function now takes the form
$$n\ln(\alpha\hat\sigma)-\sum_{i=1}^n \log\bigg[
    (1-p)\exp\bigg\{-{\big(Y_i-f({\mbox{\boldmath $x$}}_i;\Theta)\big)^2\over 2(\alpha\hat\sigma)^2}
 \bigg\}+
    {2p\over \sqrt{\pi}e}\bigg\{1+{\big(Y_i-f({\mbox{\boldmath $x$}}_i;\Theta)\big)^4
     \over (e\alpha\hat\sigma)^4}\bigg\}^{-1}\bigg]
\eqno(\number\mychapno.12)$$
The parameter~$\alpha$ can be interpreted as an estimate of how much of
the residuals in the model fit are contributed by the large
outliers generated by the contaminating distribution. If~$\alpha$
is close to~$1$ most of the contribution to the
residuals comes from the ``ordinary'' errors
associated with the normal distribution .
If~$\alpha$ is close to~$0$ most of the contribution to the residuals comes
from the large outliers associated with the fat-tailed
distribution. We expect that the maximum likelihood
estimate~$\tilde \alpha$ for~$\alpha$ will lie between~$0$ and~$1$. The parameter~$\alpha$
can also be interpreted as a ``cutoff'' switch. As~$\alpha$ is made smaller
the influence of the larger residuals on the parameter estimates 
will be reduced. This
 is much like looking at the data and removing the data points which fit
the model badly.  \XX{outliers}{in robust regression}

The code for the log-likelihood function $\number\mychapno.12$ can be
written down very concisely using \AD's vector operations. 
suppose that the observations $Y_i$ are contained in a
{\tt dvector obs} and the predicted values   $f({\mbox{\boldmath $x$}}_i;\Theta)$
are contained in a {\tt dvar\_vector pred}. The vector of
residuals $Y_i -f({\mbox{\boldmath $x$}}_i;\Theta)$ and the vector of
squared residuals can be calculated in two statments
\beginexample
  dvar_vector diff = obs-pred;     // These are the residuals
  dvar_vector diff2 = pow(diff,2); // These are the squared residuals
  v_hat = mean(diff2);
\endexample
\noindent and the estimate of the variance is given by the mean
of the vector of squared residuals.  If $\hbox{\tt b}=2p/\sqrt{\pi}e$
and $\hbox{\tt a2}=2\alpha$ the code for $\number\mychapno.12$ can be
written as 
\beginexample
  dvariable log_likelihood = 0.5*diff.size()*log(a2*v_hat)
    -sum(log((1.-pcon)*exp(-diff2/(2*a2*v_hat))
    + b/(1.+pow(diff2/(a2*v_hat),2))));
\endexample
\XX{vector operations}{example of}
\XX{vector operations}{to reduce the amount of temporary storage required}

While the combination of scalars and vectors in statements 
like {\tt 1.+v} can seem a bit strange at first one quickly 
gets used ot it. The main thing to consider is that 
some combinations of operations are not associative so that it
may be necessary to include parentheses to ensure that the 
operations are carried out in the correct order. 
As is often the case, using vector operations instead of
writing the code in terms of the components of the vector 
objects will greatly reduce the amount of temporary 
storage required.

As the examples below illustrate the parameter $\alpha$ is often not
well determined by the data. Rather than estimating this parameter
dirctly we feel that the best way to do robust regression is
probably to use  a fixed value for $\alpha$.
We recommend the value~$0.7$ for fitting the model although
we have no strong reason for preferring this particular value. 
 
\section{Using the \AD\ robust nonlinear regression routines}

The \AD\ system permits the user to obtain these robust regression
estimates in a simple fashion. Suppose that the observations $Y_i$
are contained in a {\tt dvector OBS\_Y} and that the values 
predicted by the model for these observation have been computed and
are contained in a {\tt dvar\_array PRED\_Y}. To compute the
value of the log-likelihood function $(\number\mychapno.5)$ with a fixed
value of $\alpha$ the
routine \hbox{\tt robust\_regression} is used.
The function prototype for this routine is
\beginexample
dvariable robust\_regression(dvector& OBS\_Y,dvar_vector& PRED\_Y,double& alpha);
\endexample
\noindent 
The function {\tt robust\_regression} returns minus the
value of the log-likelihood function. (The sign has been
changed for the function minimizer.)
\goodbreak
\input graph.tex
\goodbreak

To estimate the parameter $\alpha$ as well, the parameter must be passed as
a {\tt dvariable}.  The function prototype is
\beginexample
dvariable robust\_regression(dvector& OBS\_Y,dvar_vector& PRED\_Y,dvariable& alpha);
\endexample

\XX{bounding functions}{boundp}
\XX{bounding functions}{boundpin}
For the routine {\tt robust\_regression} to work properly 
when the parameter $\alpha$ is being estimated
it is necessary to
put bounds on the values that $\alpha$ can assume. This is done by
using the {\tt boundp} function supplied with \AD. The parameter $\alpha$ must
be restricted so that $0<a<=1$. (Actually the upper bound is not necessary,
but it will make the routine run a bit more efficiently.) 
The code
\beginexample
  a=boundp(x(i),.001,1.0,zpen); 
\endexample
\noindent restricts the value of $\alpha$ to lie within $0.001$ and $1.0$.
To find the initial value 
 which will give a desired initial value 
for~$\alpha$, the inverse function {\tt boundpin} can be used. To
obtain the initial value which gives a value of $0.7$ for~$\alpha$
and put the initial value in the $i$th component of the vector
of independent variables {\tt x} the code
\beginexample
  x(i)=boundpin(0.7,001,1.00); 
\endexample
\noindent would be employed.

\section{A simple example of robust regression}
\XX{robust regression} {a simple example}
We shall fit a simple 1~parameter model using
the robust regression routine {\tt robust\_regression} which
estimates the value of the parameter $\alpha$. 
The form of the model is  
$$Y_i=bZ_i+b^2$$
\noindent where the true value of~$b$ is~5.
The data for this example were produced by simulating 18 real data points
with normally distributed errors with mean $0$ and
standard deviation 40, and then 
adding two outliers which were 3.1 standard deviations and 6 standard
deviations removed from the true value. Of course such large outliers
could be identified by eye in this simple example. For large 
nonlinear models, however, there is often no simple 
graphical representation of 
the fit which can enable such a simple inspection of the residuals.
A $\circ$ was used to mark the outliers.

The results of the fitting procedure are shown in figure \number\mychapno.1.
The least squares fit with the outliers removed (estimated $b=4.45$)
is the closest to the true relationship. It represents the best estimates
which can be obtained if the outliers are identified with certainty. 
The robust fit including the outliers  (estimated $b=4.06$) did
condiderably better than the least squares fit including the
outliers (estimated $b=3.51$). 

For the robust fit the estimate of $\hat\sigma$ is $70.1$ while
the estimate of $\tilde\sigma$ is $55.6$ so that the estimate
of $\alpha$ is $0.79$. Since the real value 
of~$\sigma$ was $40.0$ the robust regression routine has underestimated
the contribution to the residuals made by the outliers. 

\goodbreak
The \AD\ code for this model is given below:
\bigskip
\listing{../../../examples/autodif/robexamp/robust.cpp}
\bigbreak
{\openup 1.3truept
\listing{../../../examples/autodif/robexamp/fcomp_r.cpp}
\bigbreak
\listing{../../../examples/autodif/robexamp/newreg2.cpp}
}
\bigbreak

\section{Setting the cutoff point in the robust regression routines}
It is important to be able to specify the ratio $\tilde\sigma/\hat\sigma$
because this parameter is often not well determined by the 
data. In the previous example the robust estimation
procedure has somewhat underestimated the contribution to the overall
variance of the residuals made by the outliers. Now 
consider the analysis for the same data as above with one change.
The value of
the second outlier has been set equal to $-75.0$ rather than $-125.0$.
Although this represents an improvement in the
quality of the data set, the parameter estimates (estimated $b=3.80$)
for the robust-mixture estimation
procedure when the parameter $\alpha$ is estimated are 
actually worse than before. The reason for this is
that the estimation of the ratio $\alpha=\tilde\sigma/\hat\sigma$ has
increased to $0.94$ which means that the robust estimation procedure
no longer recongnizes that there are outliers in the data.  
the robust-mixture estimates to the data for which the 
ratio $\tilde\sigma/\hat\sigma$
has been fixed at $0.7$ (estimated $b=4.20$) is included. 
This  is considerably better than the least squares fit (estimated $b=3.70$).
\goodbreak
\input graph1.tex
\goodbreak

We suspect that a little ``robustness'' is much better 
than none at all and that overall best performance of the 
model may be obtained by simply fixing the value for
 the ratio~$\tilde\sigma/\hat\sigma$ at some nominal value
such as $0.7$.
 

\section{A Monte Carlo evaluation of the robust regression routine}
While the above example is indicative of the superiority of the
robust-mixture estimator, it only indicates how the
estimators perform in one special case. From a statistical point of
view one is interested in the average performance of
the estimators for many different realizations of the data.
This section describes the results of such an investigation.
The data were fit by the least squares estimator, the minimum total
absolute deviation estimator and the robust-mixture estimator
(with fixed $a=0.7$). 


For each simulation a set of 30 pairs $(Y_i,z_i)$ were generated
where $z_i=i$, \hbox{for $i=1,\ldots,30$} and
$$Y_i=bz_i+b^2+\epsilon_i$$
The parameter $b$ was set equal to~5.
For the case of normally distributed model errors,
the $\epsilon_i$ are normally distributed random
variables with standard deviation 40.
For model errors exhibiting moderate deviations from normality 
the $\epsilon_i$ are a mixture normally 
distributed random
variables with standard deviation 40 with probability 0.90 and
normally distributed random variables 
with standard deviation 120 with probability 0.10.
For model errors exhibiting large deviations 
from normality $\epsilon_i$ were a mixture of normally 
distributed random
variables with standard deviation 40 with probability 0.90 and
normally distributed random variables with 
standard deviation 200 with probability 0.10.

Five hundred replications were made for each of the error distributions.
For all the simulations the value of ratio $\alpha$
was kept fixed at $0.7$ for the robust-mixture estimator. 

Now it is known that for normally distributed errors,
the usual least squares estimates are the minimum variance unbiased
estimates of the model parameters. It follows that if the
robust estimates are unbiased they must have greater variance than
the least squares estimates. This greater variance is the penalty we must
pay for using robust estimates in the case where they are not necessary.
For the problem studied here, however, the standard deviation of the
robust estimates was only about 8\% greater than the least squares
estimates for the case of normally distributed errors, while the
robust-mixture estimator had a standard deviation about 20\% smaller
than that of the least squares estimates for moderate deviations from
normality and about 35\% smaller for large deviations from normality.

\bigbreak
\vbox{\openup -2pt \tenrm
\tabskip=0ptplus100pt
\noindent Table \number\mychapno.1\quad Comparison of the least squares (LS), 
 Minimum total absolute deviation (MD), and robust-mixture (RM) estimates
for $b$ in the model \hbox{$Y=bz+b^2$.} The true value was $b=5.0$
\smallskip
\hrule
\smallskip
\hrule 
\medskip
\halign to \hsize{#\hfil\quad&#\hfil&#\hfill&#\hfill&\qquad#&
       #\hfil&#\hfil&#\hfil&\qquad#&#\hfil&#\hfil&#\hfil \cr
&\multispan3{\hfill Normal errors\quad\hfill}&&
\multispan3{\hfill Moderate deviations\quad\hfill}&&
 \multispan3{\hfill Large deviations\quad\hfill}\cr
\noalign{\kern -5truept}
&\multispan3{\hrulefill}&&
\multispan3{\hrulefill}&&
 \multispan3{\hrulefill}\cr
\noalign{\smallskip}
Mean&\quad LS&RM&MD&&\quad  LS&RM&MD&&\quad  LS&RM&MD\cr
Mean&\quad 5.01& 5.01&5.00&&\quad 4.99& 5.01&5.01&&\quad 4.98& 4.99& 5.01\cr
Std Error&\quad 0.012& 0.013&0.015&&\quad 0.016& 0.014&0.016&&\quad 0.023& 0.015& 0.016\cr
Std Dev&\quad 0.26& 0.28&0.33&&\quad 0.36& 0.30&0.36&&\quad 0.51& 0.33& 0.36\cr
\% improvement& &-0.077&-0.27&&& 20.0&0.00&& &35.3&29.4\cr
}\medskip\hrule
} %\vbox
\bigbreak
All three estimators appear to be unbiased. The standard deviation of the 
least squares estimates almost doubles as we go from the
normal errors to errors with large deviations from normality. 
The robust-mixture estimator exhibits much greater stability. 
The minimum total absolute deviation estimator performs almost as well
as the robust-mixture estimator when there are 
large outliers present,
but does not perform very well for normally distributed errors
or for medium-sized outliers present.
\bigbreak
\vbox{ \rightskip=3.2in \openup 1truept
Another way to compare the performance of the robust-mixture 
estimator over the least-squares estimator is illustrated 
by Figure~\number\mychapno.2.
The absolute value of the difference between the robust-mixture estimator 
for $b$ and the actual value of $b$ is plotted against the
absolute value of the difference between the least-squares estimator for
$b$ and the actual value of $b$. For points which lie below the
diagonal line the robust-mixture estimator has smaller absolute error
 than the least-squares estimator.
The plot is for the large-deviations set of simulations.\par}
\kern -2.85in
\hbox to \hsize{\hfill\input robgrph2.tex \hskip10pt } 
\bigbreak
\section{Investigating extremely bad data with the robust-mixture
estimator}
While the errors in the data in the previous examples exhibited moderate and
large deviations from normality, the data were not really ``bad''.
A much more difficult problem for parameter estimation is created if
all the large errors are concentrated on one side of the true
line. 
To compare the estimation methods for this kind of ``bad'' data
we generated a data set where the majority
(30 data points) of the data are distributed around the line $y=x$,
with the outliers (20 data points) distributed around the point (6,2).
The data were fit with the two parameter linear model
$$Y=a+bx$$ where the true values were $a=0.0$ and $b=1.0$.

Various fits to the data are shown in figure \number\mychapno.3.
When the value of $\alpha$ is set equal to 0.1 the estimation procedure
completely ignores the 20 outliers. We anticipate that varying
the value of $\alpha$ and noting whether the residuals suddenly shift
could be a useful way to detect these types of outliers in large
nonlinear models. Of course in a large model
it is not always possible to view the results in this simple way.
One way to examine how the vector of residuals is changing is to calculate 
the euclidean distance between the vectors. If the residuals are
changing slowly the distance between the residuals should be small.
A sudden shift in the residuals should manifest itself in a large
distance between the vectors. 
\bigbreak
$$\vbox{ \input grph3.tex } $$
\bigbreak

Table \number\mychapno.2
 shows the distance between the vectors of residuals
for different values of $\alpha$. It is clear that there are
two groups, one group for values $0.9\gt\alpha\gt0.4$ and another
group for $0.3\gt\alpha\gt0.1$. Constructing such a table could be
a possible way to search for ``interesting'' structure in the
outliers in a large nonlinear model.
\bigbreak
\vbox{
{\openup -1truept \tabskip=0pt plus 20pt \tenrm 
\noindent Table \number\mychapno.2\quad  Euclidean distance between the
vectors of residuals for the robust-mixture estimator for different 
values of $\alpha$
\medskip
\hrule
\smallskip
\hrule
\medskip
\halign to \hsize{\hfil#&&\hfil#\cr
$\alpha$&0.9&0.8&0.7&0.6&0.5&0.4&0.3&0.2&0.1\cr
\noalign{\kern -5truept}
&\multispan9{\hrulefill}\cr
\noalign{\smallskip}
0.9&0.0&\cr%&0.3&1.0&1.9&2.7&3.4&11.6&14.0&15.9\cr
0.8&0.3&0.0&\cr%0.6&1.6&2.4&3.1&11.9&14.3&16.2\cr
0.7&1.0&0.6&0.0&\cr%0.9&1.7&2.4&12.5&14.9&16.9\cr
0.6&1.9&1.6&0.9&0.0&\cr%0.7&1.5&13.4&15.8&17.7\cr
0.5&2.7&2.4&1.7&0.7&0.0&\cr%0.7&14.1&16.5&18.4\cr
0.4&3.4&3.1&2.4&1.5&0.7&0.0&\cr%14.7&17.1&19.0\cr
0.3&11.6&11.9&12.5&13.4&14.1&14.7&0.0&\cr%2.4&4.3\cr
0.2&14.0&14.3&14.9&15.8&16.5&17.1&2.4&0.0&\cr%1.9\cr
0.1&15.9&16.2&16.9&17.7&18.4&19.0&4.3&1.9&0.0\cr
%0.9&0.000&0.311&1.007&1.940&2.733&3.478&11.614&14.041&15.985\cr
%0.&0.311&0.000&0.696&1.629&2.423&3.171&11.905&14.332&16.277\cr
%0.&1.007&0.696&0.000&0.934&1.729&2.487&12.557&14.984&16.927\cr
%0.&1.940&1.629&0.934&0.000&0.799&1.574&13.424&15.848&17.790\cr
%0.&2.733&2.423&1.729&0.799&0.000&0.790&14.132&16.553&18.492\cr
%0.&3.478&3.171&2.487&1.574&0.790&0.000&14.723&17.136&19.070\cr
%0.&11.614&11.905&12.557&13.424&14.132&14.723&0.000&2.427&4.371\cr
%0.&14.041&14.332&14.984&15.848&16.553&17.136&2.427&0.000&1.944\cr
%0.&15.985&16.277&16.927&17.790&18.492&19.070&4.371&1.944&0.000\cr
} }
\smallskip\hrule} %\vbox
\bigbreak
\endchapter
\htmlnewfile 
%chapno=13
\mychapter{Problems with more than one dependent variable}

\X{several dependent variables}
\XX{several dependent variables}{use of the operator {\tt <<}}
\section{Using more than one dependent variable}
Till now all the problems considered have had only one dependent 
variable. Many problems of interest such as the problem of
solving a system of nonlinear equations require the computation
of derivatives with respect to more than one dependent variable.  

The operator {\tt <<} has been overloaded to implement 
more than one dependent variable. The user simply writes
something like
\beginexample
u << -exp(x(1))+y; // u is a dependent variable
\endexample
where {\tt u} is a {\tt dvariable} and {\tt u} will become a
dependent variable for which  derivatives can be calculated.
The operator {\tt <<} acts like {\tt =} with the additional
property that {\tt u} is made a dependent variable.
Of course it is possible to use {\tt u} in later calculations
such as
\beginexample
u << -exp(x(1))+y; // u is a dependent variable
// ...
v << u*u; // v is a dependent variable
\endexample
The process of declaring depdendent variables in this way is
global, that is, it does not have to be doen in one function.
the dependent variables are ordered in the order in which they
are encounterd when the code is executed.

The default maximum number of dedendent variables (which is at least 10)
depends on the implemntation of \AD\ which you are using.
If desired it can be increased by using the static member function
\XX{several dependent variables}{default number of}
\XX{dependent variables}{set\_NUM\_DEPENDENT\_VARIABLES(int i);}
\X{set\_NUM\_DEPENDENT\_VARIABLES(int i);}
\beginexample
gradient_structure::set_NUM_DEPENDENT_VARIABLES(int i);
\endexample
For example to set the maximum number of depdendent variables to 25
you would put the line 
\beginexample
gradient_structure::set_NUM_DEPENDENT_VARIABLES(25);
// ...  must occur before the gradient_structure declaration 
gradient_structure gs;
\endexample
\noindent into your code.
The  derivatives  are calculated by calling the function
\hbox{\tt jacobcalc(int\& nvar,dmatrix\& jac)}. The derivatives
are put into the {\tt dmatrix jac}. The minimum valid row index
of {\tt jac} must be 1 while the maximum valid row index must
be greater than or equal to the number of dependent variables.
Each row of {\tt jac} must have a minimum valid index of 1
and a maximum valid index greater than or equal to {\tt nvar}
where {\tt nvar} is the number of independent variables.


\X{gradient\_structure} \X{jacobcalc} 
\XX{finding roots of equations}{Newton--Raphson method}
\XX{Newton--Raphson method}{Roots of equations}
\XX{Newton--Raphson method}{example of}
\section {Finding roots of systems of equations with Newton-Raphson}
Suppose that you wish to solve the nonlinear system of
equations 
%$$\eqalign{
 \begin{eqnarray}
f_1(x_1,\ldots,x_n)&=0\cr
f_2(x_1,\ldots,x_n)&=0\cr
f_2(x_1,\ldots,x_n)&=0\cr
\ldots\quad&\cr
f_n(x_1,\ldots,x_n)&=0\cr
 \end{eqnarray}
%\eqno(\number\mychapno.1) $$
In general this is a very difficult problem. However, under
suitable regularity conditions, if you can
find an initial estimate for $(x_1,\ldots,x_n)$ which
is close enough then the iterative Newton-Raphson technique
will efficiently converge to the solution. Let $X$ denote the
vector $(x_1,\ldots,x_n)$
$$J(X)={\partial f_i(X)\over\partial x_j}$$
be the Jacobian matrix. Then to first order
$$ f(X+\delta)=f(X)+J(X)\delta\eqno(\number\mychapno.2)$$ where $f(X)$
denotes the vector $f_1(X),\ldots,f_n(X)$ and
$J(X)\delta$ denotes the multiplication of the vector $\delta$
by the matrix $J(X)$. Equating $(\number\mychapno.2)$ to zero and solving for $\delta$ yields
$$\delta=-J(X)^{-1}f(X)\eqno(\number\mychapno.3)$$
Now we proceed iteratively. If $X^{(r)}$ is the current estimate for
$X$ let $$X^{(r+1)}=X^{(r)} +\delta$$ be the next estimate where $\delta$
is defined by (\number\mychapno.3).
\bestbreak
\section{Implementation of the Newton--Raphson technique}
\goodbreak
This example shows the \nr\ technique for two equations. 
\bigskip
\listing{../../../examples/autodif/newt/newt.cpp}
We have chosen a rather ``degenerate'' example where the Jacobian
matrix is almost singular at the solution. The two functions are
%$$ \eqalign{
 \begin{eqnarray}
f_1(x_1,x_2)&=100x_1^2+x_2^2-1\cr
f_2(x_1,x_2)&=100x_1^2+x_2^2-1+.0001-.0001x_1^4+.01x_2^2\cr 
 \end{eqnarray}
%}$$

The calculation of the jacobian and the Newton-Raphson ``update'' step 
can be accomplished with two lines of code. 
\beginexample
jacobcalc(nvar,jacobian)
x-=inv(jacobian)*f
\endexample
\noindent The iteration continues until the norm of the {\tt dvector} {\tt f}
is less than {\tt 1.e-12}.

The only variable types occur in the
user's function {\tt function}. All the arithmetic for the
Newton-Raphson update is done with vectors and matrices which
are of constant type.  The {\tt value} function is used to
turn the {\tt dvar\_vector f} in {\tt function} into
a {\tt dvector} object so that it can be returned to the
calling routine. 

\endchapter


\htmlnewfile 
%chapno=14
\def\bs{char92}
\mychapter{A Complete Nonlinear Parameter Estimation Program}
A major difference between linear and nonlinear statistical
modelling is the inherent instability of nonlinear models.
It seems to be an empirical fact that for any
reasonably complicated model if you simply write correct
naive computer code for estimating the parameters, the
program will not work. By ``work'' we mean get the information
which you want out of the data you are trying to analyze. 
In general there are two reasons for this behaviour.
{\parindent=30truept
\begin{description}
\myitem{1.} The function being calculated is numerically unstable so
that a zero divide or a floating point overflow
 will occur, for example calculate~$1/x$ 
when~$x$ has the value~$0$ or calculate~$\exp(x)$ when~$x$ 
has the value~$90000$.

\myitem{2.} The ``correct'' values of the parameter estimates occur at a local
minimum of the objective function, not at a global minimum.
\end{description}
\par }

The problem we consider in this chapter exhibits both kinds of pathology
mentioned above.  We will describe the techniques which are necessary for
producing a stable, controllable nonlinear parameter estimation
routine for solving this problem. The code for this example is 
included on the \AD\ distribution disk.

\section{Finite Mixture Problems}
Consider a set of~$m$~random variables~$X_j$, indexed by~$j$,
with probability density functions~$F_j(x)$.
A new random variable~$Y$ can be created from the~$X_j$  as
follows. A realization of~$Y$ is produced by first picking
one of the~$X_j$ at random with probability~$p_j$ and then
taking a realization of the resulting~$X_j$. The random
variable~$Y$ is a finite mixture of the~$X_j$ with
mixing proportions $p_j$. The probability density function
for~$Y_j$ is given by~$\sum_{j=1}^m p_jF_j(x)$. 
The particular finite mixture which we shall address is
a mixture of normal distributions. Effective techniques for
estimating parameters for this problem were first considered by
Hasselblad (1966).

Assume that $X_j$ is normally distributed with mean~$\mu_j$
and standard deviation~$\sigma_j$. The probability density function
for~$X_j$ is given by
$${1\over\sqrt{2\pi}\sigma_j}\exp\Big(\,{-(x-\mu_j)^2\over 2\sigma_j^2}\Big)$$
while the probability density function for the mixture~$Y$
is given by
$$\sum_{j=1}^m {p_j\over\sqrt{2\pi}\sigma_j}
                        \exp\Big(\,{-(x-\mu_j)^2\over 2\sigma_j^2}\Big)$$
Now suppose that we are given a vector $(y_1,\ldots,y_n)$ of $n$
realizations of the random variable~$Y$. 
We wish to use the information in the~$y_i$ (and any other
auxiliary information which may be at our disposal) to estimate
the parameters~$p_j$,~$\mu_j$, and~$\sigma_j$.  

In general there are many kind of auxiliary information which
may be present in such mixture problems. For example, we may 
or may not know the
number of groups present in the mixture. We may know some or all
of the means,~$\mu_j$, or we may have some knowledge of the approximate
value of some or all of the~$\mu_j$. We may know the value of
some or all of the~$\sigma_j$, or we may know that all 
the~$\sigma_j$ have the same value.
A general routine for estimating the parameters in mixture problems
should be able, as much as possible, to incorporate auxiliary information about the parameters into the estimation process.

The maximum likelihood estimates for the parameters are found by 
finding the correct local maximum of the
log-likelihood function 
$$\myeq{\max_{p,\mu,\sigma}\,
   \bigg\{\,\,\sum_{i=1}^n\,\,\
       \log\Big[\,\sum_{j=1}^m {p_j\over\sqrt{2\pi}\sigma_j}
        \exp\Big(\,{-(y_i-\mu_j)^2\over 2\sigma_j^2}
      \,\Big)\,\Big]\,\,\bigg\}}\eqno(\number\mychapno.1)$$
\noindent subject to any constraints on the parameters which may
be appropriate for the problem. 

For mixtures of normal distributions the correct solution never occurs at 
the global maximum of the log-likelihood function when there is more
than one group in the mixture. This is due to the fact that it is
always possible to make the log-likelihood function
as large as desired by setting the mean of one of the groups equal to
the value of an observation and letting the corresponding standard
deviation parameter tend to 0. The value of the
probability density function for this group will then tend to infinity
at the point in question. It follows that to find the
correct solution the values of the standard deviations should be 
bounded away from 0.

The following code calculates the likelihood function
corresponding to equation $(\number\mychapno.1)$. Notice that 
none of the constraints which we intend to place on permissible
parameter values appears in the code. The constraint conditions
will be implemented elsewhere. Avoiding them at this stage produces
a clean piece of code which can easily be seen to
correspond to equation $(\number\mychapno.1)$. The constant terms
$1/\sqrt{2\pi}$ have been ignored. 

\goodbreak\bigbreak
\listing{../../../examples/autodif/mixture/like.cpp}

\noindent The term {\tt 1.e-20} is put in to improve the numerical stability
of the program. The routine {\tt log\_likelihood} is intended to
be used with a function minimizing routine to estimate the parameters. 
In the course of such a process the function minimizer invariably
will pick some parameters which fit certain observations extremely 
badly. An extremely bad fit is characterized by the value of
{\tt diff} being large for all groups. As a result, due to
floating point underflow,~{\tt exp(-diff*diff/(2.*v(j)))} 
can be equal to 0 for all
groups so that {\tt sum} has the value 0. Taking the logarithm
of~$0$ would cause a floating point exception or error. Adding
the small value {\tt 1.e-20} allows the program to carry on 
gracefully. 

It should be noted that Hasselblad grouped the observed data into size classes
and analyzed the grouped data. This can easily be done and should be done
if the number of observations is large, because the analysis of the grouped
data can be much faster.

\section{Putting bounds on parameter values}
\X{boundp}
\XX{bounding functions}{boundp}
\X{penalty function}
The function minimizer routine deals with a vector {\tt x}
of numbers. The likelihood calculation routine 
{\tt log\_likelihood} deals with vectors of proportions, means,
and standard deviations. We need a routine to ``put'' the x's
into the model.

\goodbreak\bigbreak
\listing{../../../examples/autodif/mixture/reset.cpp}

\X{boundp}
The function {\tt reset} not only puts the x's into the parameters of the
the function {\tt like\_fun}, it also determines which of the parameters
are active (being estimated at the time) by examining the
{\tt ivector control} and puts constraints on the values the
parameters through the function {\tt boundp}.
The function
\beginexample
dvariable boundp(dvariable x,double fmin,double fmax,dvariable& fpen)
\endexample
\X{penalty function}
\noindent bounds the variable {\tt x} to the closed interval 
$[${\tt fmin,fmax}$]$.
The term {\tt fpen} is a penalty term which modfies the behaviour of the
function. It must be added to the function being minimized.
{\bf Note that this means that if you want to maximize a function you must
change the sign of the function before adding {\tt fpen} to it!} 
The function {\tt reset} supplies complete control over which parameters
are active through the vector {\tt control}.  For example,
if {\tt control(1)} is equal to 1 while {\tt control(2)} and 
{\tt control(3)} are equal to 0, then only the proportions at age will
be estimated while the means and standard deviations will be kept fixed.
The {\tt dvectors mumin} and {\tt mumax} determine lower and upper
bounds on the means while 
the {\tt dvectors sdmin} and {\tt sdmax} determine lower and upper
bounds on the standard deviations. 

Suppose that you want the standard deviations to be the same for all the
groups. This can be accomplished easily by modifying the last part of 
{\tt reset} to
\beginexample
  if (control(3)>0)  // control(3) determines whether the sd's are active
  {
    for (int j=1;j<=ngroups;j++)    
    {
      sd(j)=boundp(x(ii),sdmin(j),sdmax(j),zpen); //Note that ii is not
                                                  // incremented
    }
    ii++;   // Now increment ii
  }
\endexample
Another modification would be to enable individual mean or standard 
deviations to be fixed or active.
\section{Getting the initial x vector}
\X{boundpin}
\XX{bounding functions}{boundpin}
It should be clear that the vector {\tt x} is a transitory object
associated with the minimization routine. Without the function {\tt reset}
it is not even clear to which model parameter a particular component
of {\tt x} corresponds or what value of the parameter it produces.
For this reason it does not makes sense to save the values of
{\tt x}  during or after the analysis. It is the model
parameters {\tt p}, {\tt mu}, and {\tt sd} which should be saved.
However this raises the question: Given the values of the
model parameters, the values of the control vector, and the
values of the bounds on the parameters, what should the initial values 
of the {\tt x} vector be? This problem is solved by the functions
{\tt boundpin} and {\tt xinit}.

\bigbreak
\listing{../../../examples/autodif/mixture/xinit.cpp}

The control structure of {\tt xinit} exactly parallels that of {\tt reset}
so that no matter what switch combinations or bounds are used, the
correct values will be put into the components of {\tt x}.
The function
\beginexample
double boundpin(double x,double fmin,double fmax)
\endexample
\noindent is the inverse of the function {\tt boundp}.
Notice that the control vector, {\tt control}, determines how many 
active parameters there are. 

The job of the function
{\tt nvarcal} is to use this information to calculate the number
of active parameters.
\bigbreak
\listing{../../../examples/autodif/mixture/nvarcal.cpp}

The three functions {\tt reset}, {\tt xinit}, and {\tt nvarcal}
all share a common control structure and make it possible to
coordinate the various users options in the estimation routine.

Here is the code for the main routine in the mixture analysis 
program. The data are assumed to be in a file named {\tt mixture.dat}.
The initial values for the parameters are assumed to be in a file
named {\tt mixture.par}.
\bigbreak
\listing{../../../examples/autodif/mixture/mixture.cpp}

The job of the function fcomp is to put the x values into the model
by calling the function {\tt reset} and then to evaluate the log-likelihood
by calling the function {\tt log\_likelihood}.
\bigbreak
\listing{../../../examples/autodif/mixture/fcomp_m.cpp}

Notice that while the parameterization of the mixture proportions,~{\tt p},
has restricted the proportions to lie between~$0$ and~$1$, the
proportions have not been restricted so that they sum to 1. 
The function {\tt normalize\_p} restricts the~{\tt p} so that they sum 
to~$1$. 
\bigbreak
\listing{../../../examples/autodif/mixture/normaliz.cpp}
\X{penalty function}
\noindent The penalty {\tt 1000.*log(psum)*log(psum)} has been added to the
function being minimized to remove the degeneracy in the parameterization of
the~{\tt p}. Without this penalty the minimizing value 
would be independent of the value of {\tt psum}. Including
the penalty ensures that the minimum value will occur for the value
{\tt psum}$=1$.
\section{Saving the parameter estimates}
\AD\ uses the same memory locations for the derivative
calculations as are used for holding the variable objects.
In consequence, after the derivatives have been calculated
the current value of the parameter estimates has been lost.
In order to save the results after the minimization has
been carried out it is necessary to calculate the
parameters one more time and print them into a file.
The simplest way to do this is to use the code for
{\tt fcomp} which already exists, and to put a
switch into {\tt fcomp} so that saving the parameters is enabled.
The parameters are written into the file {\tt mixture.par} by the
routine {\tt save\_pars}.
\bigbreak
\listing{../../../examples/autodif/mixture/savepar.cpp}
\endchapter


\htmlnewfile 
%chapno=15
\mychapter{A Neural Network}
\X{neural net}
\XX{neural net}{feed-forward}
This example illustrates the use of \AD\thinspace's two and three dimensional
ragged arrays.

\section{Description of the neural network}
An $n$-layer feed-forward neural network can be viewed as a process which, 
given a set of $m_1$ inputs $(x_{11},\ldots,x_{1m_1})$ at layer $1$, 
produces a set
of $m_n$ outputs $(x_{n1},\ldots,x_{nm_n})$  at layer $n$.
Level $1$ is referred to as the input layer. Level $n$ is
referred to as the output layer. Levels $2$ through $n-1$ are
referred to as hidden layers. The quantities $x_{ij}$ are referred to
as nodes.
\XX{neural net}{input layer} \XX{neural net}{output layer} \XX{neural net}{nodes}

In a feed forward neural network the values of the nodes in each
layer are determined by the values in the previous layer via the
relationship
$$\myeq{x_{i+1,j}=\Phi\Big(\sum_{k=1}^{m_{i-1}}w_{ijk}x_{ik}+b_{i+1,j}\Big)}
  \eqno(\number\mychapno.1)$$
The function $\Phi$ is referred to as a squashing function.
To the best of our knowledge, the particular
form of $\Phi$ does not seem to be especially
important so long as it is a continuously differentiable monotone 
increasing function
which maps $(-\infty,\infty)$ into some desired bounded set $(a,b)$. 
The bounded set is often of the form $(-0.5,0.5)$ or $(0,1)$.
The $w_{ijk}$ are referred to as weights, while the $b_{ij}$
are called bias terms. 
\XX{neural net}{squashing function} \XX{neural net}{bias terms}

The shape or topology of the neural net is completely determined
by the number of levels, $n$, and the number of nodes in each
level, $m_i$. Assume that the $m_i$ are contained in an 
{\tt ivector} object {\tt num\_nodes}.

\section{Implementation of the neural network}
\X{two dimensional ragged array}
The $x_{ij}$ can be viewed as a two dimensional ragged array
with $n$~rows and $m_i$~elements in the $i\,$'th~row. Such
an array can be created by the declaration
\X{ragged matrix}
\beginexample
dmatrix x(1,num_levels,1,num_nodes); // these are the nodes
\endexample
\noindent The bias terms $b_{ij}$ can also be viewed as
a ragged matrix {\tt B} with $n-1$ rows. The declaration for B
requires a slight modification of the {\tt ivector num\_nodes}
to reflect the fact that {\tt B} does not need to contain
row 1. The legal row index values for {\tt B} should be from 2 to $n$.
\beginexample
ivector iv(2,num_levels);
for (int i=2;i<=num_levels;i++)
{
  iv(i)=num_nodes(i); // Put the desired elements of num_nodes into iv
}
// Now declare B
dmatrix B(2,num_levels,1,iv);
\endexample
For each $i$, $1\le i<n$, the weights, $w_{ijk}$,  can be viewed as
an $m_{i+1}$~by~$m_i$~matrix.
Thus the $w_{ijk}$ can be viewed as a ragged three dimensional array {\tt W}.
To create {\tt W} we need some
{\tt ivector} objects to contain the variable number of rows and columns.
\X{three dimensional ragged array}
\beginexample
ivector nrows(1,num_levels-1);
ivector ncols(1,num_levels-1);
for (i=1;i<num_levels;i++)
{
  nrows(i)=num_nodes(i+1);
  ncols(i)=num_nodes(i);
}
// Now declare W
d3_array W(1,num_levels-1,1,nrows,1,ncols);
\endexample
If the squashing function $\Phi$ is
extended to vector objects by acting elementwise on the elements
of the vector, the code segment to carry out the calculations in 
equation $(\number\mychapno.1)$
can be written as
\beginexample
for (int i=1;i<num_levels; i++)
{
  x[i+1]=Phi(W[i]*x[i]+B[i+1]);
}
\endexample
\section{ Initializing the vector of active parameters and resetting
active parameter values with the set\_value functions}

To set up the initial vector of parameter estimates  needed by the
function minimizer it is necessary to copy all the values of the active
parameters. To copy values into or from two and three dimensional objects
one must use a series of nested loops. This process is tedious and prone
to error. A collection of functions has been included in \AD\ to
carry out these sort of tasks. The function 
\beginexample
set_value_inv(A,y,ii); // this is the ``inverse'' of the set_value 
                       // function
\endexample 
\noindent will copy the elements of {\tt A} into the {\tt dvector y}
starting at the offset specified by the {\tt int ii}. {\tt ii}
will be incremented so that it points at the
correct offset in {\tt y} to add the values of other objects.
{\tt A} may be a {\tt dvector}, a {\tt dmatrix}, or a
{\tt d3\_array}.
For a {\tt dvector v} the statements
\beginexample
int ii=6;
set_value_inv(v,y,ii);
\endexample
\noindent will have the same effect as
\beginexample
int ii=6;
for (int i=v.indexmin();i<=v.indexmax();i++)
{
  y(ii++)=v(i);
}
\endexample
\X{set_value_inv}\X{set_value}\X{size_count}
\X{Counting the number of elements in container classes}
The function
\beginexample
set_value(A,y,ii);
\endexample
\noindent takes the values from the {\tt dvar\_vector y} beginning
at the offset {\tt ii} and
puts them into the object {\tt A}.
{\tt A} may be a {\tt dvar\_vector}, a {\tt dvar\_dmatrix}, or a
{\tt dvar3\_array}. Of course it is the users responsibility to declare the
vector objects {\tt y} to have the appropriate size.
To count up the number of elements in the container objects the
function {\tt size\_count()} has been defined.
The functions {\tt set\_value()}, {\tt set\_value\_inv()}, and 
{\tt size\_count()} have been used in the neural net code 
to simplify the construction of the functions 
{\tt nvar\_calc} which calculates the
number of active parameters to be estimated, {\tt yinit} which puts the
intial value of all the active parameters into the vector {\tt y} for
the function minimizer,
and {\tt reset} which calculates 
the current estimate of the active 
parameters from the vector {\tt y}. The code for these functions follows.
\beginexample
// count up the number of active parameters
int nvar_calc(dmatrix& B, d3_array& W)
{
  int nvar=0;
  nvar+=size_count(B);
  nvar+=size_count(W);
  return nvar;
}

// Put the parameter values from W and B into the vector y
void yinit(d3_array& W,dmatrix& B,independent_variables& y,int& num_levels,
  ivector& num_nodes)
{
  int ii=1;
  set_value_inv(W,y,ii);
  set_value_inv(B,y,ii);
}

// Put the values from the dvar_vector y into W and B
void reset(dvar3_array W,dvar_matrix& B,dvar_vector& y,int& num_levels,
  ivector& num_nodes)
{
  int ii=1;
  set_value(W,y,ii);
  set_value(B,y,ii);
}
\endexample

\section{Training the Network}
\XX{neural net}{training set} \XX{neural net}{weights} 
\XX{neural net}{bias terms}
A training set for a neural net is a collection of input layer vectors
together with a corresponding set of output layer vectors.
For each input layer vector the corresponding output layer
vector can be thought of as the desired correct output response of the 
network for the given input stimulus. The weights and bias terms of the
network are adjusted so that the difference between desired response 
and the actual response is minimized.
We shall assume that the training set is contained in a file
called {\tt learning.smp}. The training set  will be read in every time the
response of the neural net is calculated. This provides
for very large training sets which could not be kept in memory.
\beginexample
ifstream infile("learning.smp");  // Open the file with the training data
infile >> num_learn;  // read in the number of examples to be used
// Loop through the training data 
for (int k=1;k<= num_learn; k++)
{
  infile >> x[1];
  if (!infile)
  {
    cerr << "Error reading x[1] from file learning.smp\n";
    exit(1);
  }
  infile >> correct_response;
  if (!infile)
  {
    cerr << "Error reading correct_response from file learning.smp\n";
    exit(1);
  }
  for (int i=1;i<num_levels; i++)
  {
    x[i+1]=Phi(W[i]*x[i]+B[i+1]);
  }
  dvar_vector vdiff=correct_response-x[num_levels];
  z+=vdiff*vdiff; // z sums up the squared difference between the
                  // desired response and the actual response

  // Add the penalty terms
  dvariable v;
  v=norm(W);
  z+=.01*v*v;
  v=norm(B);
  z+=.01*v*v;
}
\endexample
The parameters for {\tt W} and {\tt B} are chosen so that the
total squared error is minimized. A small quadratic penalty
function of the model parameters is added to keep the values of
the parameters close to 0. This is much like the use of ridge regression
(Hoerl and Kennard, 1970)
in linear models.  We don't know whether the use of such penalty
functions really improves the performance of the neural 
net, but it seemed to produce good results for the
 problem considered here.
\X{ridge regression}
\section{The squashing function}
For the squashing function, $\Phi$, we have employed the inverse tangent 
function~{\tt atan}.
\XX{squashing function}{use of inverse tangent for}
\beginexample 
dvar_vector Phi(dvar_vector& v)
{
  // This squashing function maps into the interval (-.5,.5)
  RETURN_ARRAYS_INCREMENT;  // Need this because the function returns a 
                            // variable type
  int min=v.minindex();
  int max=v.maxindex();
  dvar_vector tmp(min,max)
  for (int i=min;i<=max;i++)
  {
    tmp(i)=0.31831*atan(v(i)); //0.31831 is 1/PI 
    // tmp(i)=0.31831*atan(v(i))+0.5; // This maps into (0,1) 
  }
  RETURN_ARRAYS_DECREMENT; 
  return(tmp);
}
\endexample
\section{Specifying the structure of the neural net and getting initial parameter estimates}

The structure of the neural net is determined by specifying the number of
layers and the number of nodes in each layer. This information
is put into the file {\tt struct.nrl}.

To begin training the neural net the {\tt d3\_array W} and the
{\tt dmatrix B} were filled with normally distributed random numbers,
with mean 0 and standard deviation 1.0.  The value of 1.0 for the standard
deviation was picked by trial and error. If the standard deviation is
too small the minimization routine doesn't seem to be able to find 
a good minimum and the procedure ``stalls'' at a bad local minimum.

The program {\tt WINIT.EXE} reads the file {\tt struct.nrl} and
writes the initial parameter values into the file {\tt WEIGHTS.PAR}.  
 

\section{A three spiral problem}
Touretsky and Pomerlau (1989) discussed the problem of training
a neural network to recognize points which lie on or  near
a double spiral. We thought it would be interesting to
extend this problem to points which lie on or near 
a triple spiral.
\bestbreak
\vbox{ \rightskip=2.7in \tolerance=5000 \openup-.5truept
The training set for the neural net is shown in Figure \number\mychapno.1.
The net is being trained to determine whether a point in a disk 
lies on or near one of three spirals. 
The input level to the net consists of two nodes, the 
cartesian coordinates of the point. The output level of the net has 
three nodes. The node corresponding to the spiral on which
a particular point is located
is intended to be turned on (have value +.5).
 and to be turned off (have value -.5) otherwise. Actually,
any value for a node that is less than 0 is considered
``off'' while any value which is greater than $0$ is considered
``on''. Due to the nature of the inverse tangent squashing function
the maximum value produced by the neural net is about $0.42$ while the
minimum value is about $-0.42$. If the neural net has $0$, $2$, or $3$ output
nodes turned on ($>0$) for a point then that point is considered unclassified.
\par} 
\kern -2.8in
\hbox to \hsize{\hfill\input grph5.tex \hskip10pt } 
\bigbreak
\vskip \baselineskip
A typical input consist the pair of coordinates of a point
on the spiral $(0.0618,0.1901)$ together
with the triple $(0.5,-0.5,-0.5)$  which is the desired output for 
this point since it lies on the first spiral arm.

We trained neural nets of increasing complexity until a design 
was obtained which had enough parameters to classify all the points
in the learning sets correctly. The neural net we chose
has 6 layers with 2 input nodes, 8 nodes in each of the 4 hidden layers,
and 3 output nodes. This neural net has 267 parameters to be estimated.

\X{GRADFIL1.TMP} 
Training the neural net required about 20 hours on a 20 MHz 386
computer. The file {\tt GRADFIL1.TMP} which stores the temporary
information needed to calculate the gradient was about 3 MB in size.
It was stored on a RAM disk. At any time it is possible to stop
the program by typing ``q''. The current parameter estimates
are printed into the file {\tt EST.PAR}. To restart the learning
process you must first copy the file {\tt EST.PAR} into the
file {\tt WEIGHTS.PAR}. In this manner it is possible to 
carry out training sessions during periods when the computer is not required.
\bigbreak
\input grph6.tex
\bigbreak
Figure \number\mychapno.2 illustrates how the neural net has classified the 
points in the disk. The points in A, B, and C have been classified
as being associated with the corresponding spiral. The small points in D
are unclassified. The large points in D are the training set. The unclassified 
points appear to mark fairly regular boundaries between the spirals. 
Our solution does not seem to exhibit the ``bumps and gaps''
which were reported by Touretsky and Pomerlau for the double spiral solution.
Perhaps this is due to the superiority of our minimization
procedure over the method they employed.
It may also be due to the use of more parameters combined with a 
small quadratic penalty 
function on the parameters to introduce stability into the model.
In any event we seem to have obtained a nice solution
to the problem.
\vfill\eject
\input auto_doc.bib
\endchapter
\pagestyle{xxx}
\input all_au.ind
\end{document}
